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NUMERICAL APPROXIMATION 
The book gives a careful discussion of selected methods 
for such numerical operations as tabulation, interpolation, 
differentiation, integration (quadrature), and curve fitting. 
Each process is illustrated by an actual numerical example 
complete with working detail and with notes on features 
of importance. 

The selected numerical processes are embedded in a 
discussion of the approximate representation of tables of 
numbers and of ‘complicated’ functions over limited 
intervals using polynomials of low degree. Such repre¬ 
sentations form a natural basis for the development by 
finite-difference and least-square techniques of the various 
numerical processes which are introduced. 

An approach of tlris kind, in which particular numerical 
processes are set against a proper mathematical back¬ 
ground, provides students with a reasonable perspective of 
numerical analysis as a branch of mathematics and not a 
mere catalogue of formulae. 

‘Anyone who is concerned with numerical computation 
or who is interested in the application of mathematics in 
such a way as to extract from a table of figures the most 
precise information that it is capable of giving will find 
this book useful.’ — Engineer 
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Preface 


Numsucal computation plays an important and growing 
role in the physical sciences, engineering and technology! 
and it serves across the whole range from the development 
of new theories to the analysis of experimental data. As 
problems increase in complexity and as high-speed com¬ 
puting machines spread so that their use is the rule rather 
than the exception, it becomes increasingly urgent that 
most students in these subjects should have some introduc¬ 
tion to numerical techniques* 

The effective use of computing machines involves two 
stages* First we need an understanding of the kinds of way 
in which numerical processes can be constructed and put to 
use; and then we must learn the system of programming 
for the particular machine at hand. Programming systems 
are like languages which enable u$ to describe problems in 
terms acceptable to the machine; at best they are remark¬ 
ably easy to learn, but In the long run their effective use 
must depend on a proper preparation of problems and 
hence on a thorough understanding of numerical analysis. 

In this short book (and a forthcoming companion 
volume Numerical Solution of Equations) I have chosen to 
concentrate my attention on the mathematical aspects, be¬ 
cause my experience has shown that in the end students 
will master the use of fast machines more quickly if the 
early stages of teaching are separated, so that they are intro¬ 
duced first to the ideas underlying numerical processes, and 
then to their adaptation to machine use. 

Like any other branch of mathematics, numerical 
analysis must be understood if it is to be of any real value, 

V jj 
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For this reason I have discussed relatively few of the many 
formulae available (as long lists of formulae add nothing 
to understanding), and I have avoided the use of opera¬ 
tional methods (as these often confuse inexperienced 
students by obscuring the significance of steps in the deriva¬ 
tion of a formula). 

The examples and exercises are particularly important, 
as it is difficult fully to grasp numerical processes in the 
abstract; they form a real part of the book, and the reader 
is recommended to follow all numerical steps using such 
calculating aids as may be available. Nothing is simpler 
than to construct fresh exercises from any set of tables. 

B. R. MORTON 

The University, 

Manchester 


CHAPTER ONE 

Introduction to Numerical Analysis 


1 , 1 . The need for numerical methods of analysis 

Progress in the physical, engineering and other sciences 
is based on two parallel but strongly coupled streams. On 
the one hand experimental, involving the design of experi¬ 
ments and the collection of physical data; and on the other 
hand theoretical, through the construction of model theories 
and their use in the interpretation of experimental results 
and in planning new experiments. The results of experi¬ 
ments are obtained as numbers or as tables of numerical 
values, and most theoretical solutions have also to be ex¬ 
pressed numerically before they can be assessed critically 
by comparison with experiment. 

The most important mathematical techniques used in the 
development of model theories are the normal processes of 
analysis, many of which are introduced in other books of 
this series. These are powerful processes which are of the 
greatest possible importance to the applied mathematician, 
but they often produce formal solutions rather than the 
final numerical solutions that are needed for immediate 
display of the physical significance of a solution and for 
direct experimental comparison. For example, it is easy to 
show that the solution to the differential equation 

dyjdx+ 2 xy= 1 with y = 0 at *=0 is y =e - * 2 e l 'dt, 

but this formal solution is of little numerical value as it is 
easier to solve the original equation by a direct numerical pro¬ 
cess than to carry out a numerical integration of the rapidly 
growing integrand (The ‘formal’ solution is numerically 
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useful given tables of e~ x * and 



e t 2 di t but these will have 


been calculated numerically and the use of such tables is 
itself a numerical process; moreover, this particular in¬ 
tegral is quite likely to have been tabulated through a direct 
numerical solution of the differential equation!) There are 
other types of problem, such as the solution of non-linear 
differential equations, which are of central importance in 
current work on physical theories but which cannot yet be 
handled analytically to any general extent, and numerical 
methods of solution have an essential role to play here if 
theory is not to lag seriously behind experiment. Again, the 
normal methods of analysis cannot be applied directly to 
operations on tables of experimental values or other tables, 
and such operations must obviously be carried out using 
numerical processes* 

Numerical work will often enter as a normal part of the 
mathematical treatment of a problem, and numerical 
methods form an essential part of the standard equipment 
of any applied mathematician* 


1-2- Approximations to functions 

A basic analytic method is the expansion of a 'compli¬ 
cated 1 function as a convergent infinite series of terms with 
simpler (or otherwise more appropriate) functional form. 
These expansions may be regarded as infinite approxima¬ 
tion processes in which the error may be made arbitrarily 
small by taking progressively more terms into account. 

The most important case is that of Taylor expansion in 
an infinite power series f where the function f{x) is repre¬ 
sented for some neighbourhood of a chosen point xq by the 
infinite series, 

f( x ) =/( A 'o) + (# — xo)f'(xo) + - Wo) 2 /"(.ro) + 

... + ^j(ac - x 0 ) n fW(xo) + ... 
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APPROXIMATIONS TO FUNCTIONS 


in powers of the displacement x — where /' — dfjdx, 

f(n) = d n fjdx n y and/<*>(*&) is to be calculated at x-xq. Such 
an expansion will be valid within some range of con¬ 
vergence tf(xo) -Oc <6(^o), and within this range the trun¬ 
cation error produced by discarding all terms after the rath 
can be made smaller in magnitude than any given positive 
constant by judicious choice of n. The difficulty from a 
numerical point of view is that this convergence may be so 
weak as to be useless in practice, and also it describes the 
ultimate behaviour of the terms in the series whereas 
numerically we want to concern ourselves with only the 
leading terms* Thus, the expansion 

tan -1 x=x — —+ 


is convergent for — 1 ^ x ^ 1; numerically it is very effec¬ 
tive for small values of x and only two terms are needed at 
£_{H for five decimal accuracy, but the number of terms 
needed increases rapidly with #, to seven terms at jc — 0*5 
for five decimal accuracy and fifty-one terms at x=l for 
two decimal accuracy* 

The Taylor series can also be written in the terminated 
form 

/(*) =f(x o) +(* - Xo)f'(xo) + •••+£](*-+*»• 


with remainder R n - 


1 

0»+i)i 


—s&)* +1 / (n+1 >( £) for some f in 


the interval {x 0 , #). The remainder R n may be regarded as 
the error involved in the approximate representation 


1 

f(x) —/(#o) + (* - xo)f(xo) + ... + “[(*- #o) n / <B) (*°)> 


and it may be seen that this is a simple polynomial repre¬ 
sentation of degree n for the function f{x). But can we 
reasonably expect that such polynomial representations can 

3 








INTRODUCTION TO NUMERICAL ANALYSIS 

in practice be used without substantial error caused by 
neglect of the remainder R n ? This question is answered 
partially by Weierstrass’s theorem on approximation by 
polynomials. 

Weierstrass’s theorem* states that a function continuous 
in a closed interval can be approximated over the whole 
interval by a polynomial of suitable degree which differs from 
the function by less than any given positive quantity e at 
every point of the interval . As the allowable error of approxi¬ 
mation e is decreased, the degree of the approximating 
polynomial mast be increased. This theorem is of the 
greatest importance because it implies that in practice, 
where we work always to a limited and specified accuracy 
only, we can replace any continuous function over a given 
range by an approximating polynomial without significant 
loss of accuracy* Hence we are justified in assuming a poly¬ 
nomial approximation for a given set of experimental or 
other values, and also in matching given functions by poly¬ 
nomial approximations over closed intervals* There re¬ 
mains one major drawback, in that Weierstrass f s theorem 
gives no indication whatever of the degree of polynomial 
required for a given accuracy of approximation. This is 
typical of the difference in attitude in pure and in numerical 
analysis; the pure analyst is concerned first with the 
general principles of approximation, but the numerical 
analyst wants to know when he can reasonably use quad¬ 
ratic or cubic approximation! We shall see that in practice 
there are ways of judging the degree of a suitable poly¬ 
nomial representation even when no direct estimate for the 
remainder l?n can be made* 

The way is now open for a major change in emphasis. 
We shall adopt unashamedly the idea of simple polynomial 
approximation to the functions or tabular values involved 
in calculations as a basis for developing numerical processes. 
Moreover, we shall reinterpret the significance of the re- 

* For a proof and various extensions see e.g* Jeffrey* and Jeffrey* 
Methods of Mathematical Physics, Cambridge, 
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1.2 APPROXIMATIONS TO FUNCTIONS 

mainder or error term of the terminated Taylor expansion 
as providing a guide to the working range over which a 
polynomial representation of specified degree can be used to 
the accuracy required. We shall thus construct numerical 
processes in which each function is represented by a poly¬ 
nomial of reasonably low degree (so that the ‘unit opera¬ 
tions’ may be easily carried out) and in which the range of 
calculation is subdivided into working intervals which are 
short enough locally to preserve the desired accuracy over 
the whole calculation. 

There are also many types of infinite series expansion in 
terms of complete sets of orthogonal functions. The most 
familiar of these are Fourier series* in terms of cosines, 
sines, or a mixture of both, and in each case valid within a 
range of orthogonality of the set of functions involved. 
Other sets of orthogonal functions which may be used for 
series expansion in an analogous manner include Legendre 
polynomials, Chebychev polynomials, Bessel functions, and 
so on. In principle, each of these expansions can be used in 
truncated form as a basis for the development of numerical 
processes in very much the same way as direct or 'simple 1 
polynomial representation. In practice, some of the most 
useful methods have been developed from truncated 
Fourier series and series of Chebychev polynomials. 
Fourier series methods form the basis of harmonic analysis 
and are particularly useful in handling functions or data 
with natural periodicities. Chebychev polynomial methods 
provide the most rapid convergence of all approximations 
based on special or simple polynomials, and give a special 
degree of uniformity in the distribution of errors of ap¬ 
proximation over the range of application (in contrast with 
simple polynomials where the error of approximation tends 
to be smaller in the centre and to increase fairly sharply 
towards the ends of the range). It may be noted that a 
similar form of Weierstrass’s theorem holds for Fourier and 
Chebychev representations. 

See Sneddon, Fourier Series , in this serie* 
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All the methods of approximation which have been 
described can be collected into the single expression 

f(x) o(a?) + + ... + OnM*)> 

and any such finite approximation may be used as a basis 
for the development of numerical processes (though only 
with complete sets of functions can we be sure that the error 
will be reduced progressively as we include successively 
more terms of the representation)* This representation is an 
approximate form with error /?„(*) of the exact terminated 
expansion 

f(x) = ao<f>o(x) + tf> i(x) + ... + a n <f>n(x) + R»(x). 

In order to emphasize the role of approximation in the 
processes we shall develop, the symbol will be used to 
denote any of these approximate representations, such as 
polynomial approximation with <f>i(x) = x k or Fourier 

approximation with ^*(#) = s i n Many authors use the 

symbol = for all purposes, but this can prove confusing to 
the beginner. In each case the maximum value of j? n (a:) (if 
it can be estimated) in a range of application provides a 
measure of the accuracy obtainable with that approxima¬ 
tion. 

A family of numerical processes will be developed in the 
following chapters from the simple basis of polynomial 
approximation, using the general form 

f(x) ^ +aix +*2*® + ... + a n x n , 

(with appropriately chosen degree n) for known and un¬ 
known functions and for sections of tabulated values. 
Numerical processes based on polynomial representation 
are relatively easy to derive, simple to apply, and very satis¬ 
factory for displaying the principles of numerical approxi¬ 
mation ; moreover, they require only a very modest mathe¬ 
matical background, and so form an ideal introduction to 
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I .3 NUMBERS AND TABLES OF VALUES 

the ideas underlying numerical analysis. In contrast, the 
development of methods based on approximation by special 
functions calls for a wider knowledge of analysis; hence in 
spite of their greater current importance, methods based 
on approximation by orthogonal functions will be post¬ 
poned to a later book in the series, 

1.3. Numbers, continuous functions, tables of values, 
and tabulated experimental results 

These provide the raw material of numerical analysis. 

We shall normally represent numbers in decimal form.* 
Thus the number 

3*14159 

has five decimals (or decimal places, or decimal digits), and 
six significant figures. It is the value of tt, and is to be 
interpreted as being correct to within five units in the sixth 
decimal place; it is impossible to guarantee greater accuracy 
without writing additional decimals, and it is an important 
convention that it should be regarded as dishonest to fill 
decimal places with digits known to be inaccurate (or, 
indeed, with digits of unknown accuracy). This represents an 
approximate value of v which is of course an irrational 
number, but irrational numbers cannot as such enter 
numerical analysis since we are forced to approximate to 
them. In the same way we must approximate all rational 
numbers which cannot be expressed exactly in the number 
of decimals with which we choose to work, and for obvious 
practical reasons this number cannot normally be very 
large. 

Sometimes we have to handle very large or very small 
numbers, such as Avogadro’s number 6 023 x IQ 23 , or the 
charge of the electron in e.s.u. 4-802 x IQ -10 ; and in such 
cases we usually separate the significant figures (which we 
write in the form of a workable decimal number of order 


# This U by no mean® the only possibility, though it js most con¬ 
venient for our purposes; modern computing machine® uae binary form. 
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unity) from the powers of ten involved, and we carry the 
two parts of the number separately through the calculation. 
This may be called floating point decimal form, and may also 
be written 6 - 023 , 23 and 4 - 802 , —10 for these examples. 
It may be noted that the value for the charge of the electron 
implies that it is known correctly to about one part in ten 
thousand (to be precise, 5 parts in 4S020). 

Continuous functions could in principle be represented 
numerically by curves showing their value at every point, 
but for practical purposes the only way in which most 
functions can be represented accurately and conveniently is 
by a table of values at specified intervals of the indepen¬ 
dent variable,* And of course we are only too familiar with 
such tables of logarithms, sines, cosines, exponentials and 
so on. In these printed tables each entry represents a true 
value, correct to within half a unit in the last digit unless 
someone has blundered, of a precisely defined function. (We 
must accept the fact that other people make mistakes just 
as we do ourselves, and quite a lot of mistakes appear in 
print!) These functions obviously exist between tabulated 
values, and we shall develop the process of interpolation to 
calculate such intermediate values. However, we cannot 
tabulate the functions more closely and many functions are 
tabulated scarcely at all because of the sheer magnitude of 
the task involved. Thus, a table of five-figure logarithms 
for the numbers from 10 by intervals of 0-1 to 100 occupies 
two pages with 900 main entries; a similar table of a func¬ 
tion of two variables would take two large books of 900 
pages each; one of three variables a room of books; one of 
four variables (such as the hypergeometric function) a 
library of more than one and a half million volumes! The 
evaluation of many functions can be carried out very 
quickly on fast computing machines. 

Tables of experimental results differ from tables of functions 

* But see Jabnke and Emde, Tables of Functions, Dover for some 
beautifully drawn curves which are extremely effective in displaying the 
broad types of behaviour of rather complicated functions. 
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in that each entry is likely to include experimental error, and 
although we can place probable bounds on the magnitude of 
the error we can do nothing to improve the accuracy of indi¬ 
vidual results without repeating and refining the experi¬ 
mental measurements. There are two kinds of experimental 
error: systematic and random. Systematic errors may be due 
to a persistent fault in the performance or faulty design of 
the experiment; they can seldom be detected by direct 
analysis of the experimental results as each measurement 
suffers, and they are less a part of the analysis of the results 
than of experimental design. Random errors arc often ap¬ 
parent on inspection, specially when the results are pre¬ 
sented graphically and show ‘scatter’; they can be kept at a 
low level by careful work but they cannot be eliminated, 
and they set an obvious limit to the accuracy which can be 
obtained in calculations using the experimental results. 

A common way of approximating a set of experimental 
results is to plot them on suitable logarithmic scale or other 
graph paper in an attempt to produce a linear relationship; 
a ‘line of best fit’ can then be estimated by eye, and will 
serve to define a functional approximation to the original 
results. When the scatter is appreciable, there is some temp¬ 
tation to use one of the smoothing techniques so that the 
adjusted points lie more closely on a smooth curve; how¬ 
ever, such devices should be used with great caution, as the 
raw experimental results provide our only direct link with 
the physical world and if they do not form an obvious 
pattern this may be because we have taken the wrong 
measurements. We shall mention more sophisticated 
methods of curve fitting in chapter 5 , and it should be noted 
that there are statistical methods for assessing the signifi¬ 
cance of variations in experimental measurements. 

1.4. Errors and accuracy in numerical analysis 

The only useful numerical solution is a correct one— 
though we must be careful how we interpret the word 
correct’, We can usefully distinguish two causes producing 

8 9 
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wrong numerical answers: errors which arise from the finite 
character of numerical analysis and are inherent in our 
processes, and blunders caused by our own carelessness of 
the failure of machine aids. 

Inevitably we must work with a limited number of sig¬ 
nificant figures, because our machines have limited capacity 
and indeed we have even more limited capacity ourselves 
in practice 1 Few numbers will be small integers, and there 
will usually be a rounding off error produced by the rejection 
of all digits after a selected decimal place. Our procedure 
for rounding off will be to preserve unchanged the last 
digit retained on discarding fewer than five units in the next 
place, to increase it by one on discarding more than five 
units, and to round off to the nearest even number on dis¬ 
carding exactly five units. Some authors always round five 
units in one direction (usually up), but our procedure 
preserves a randomness and so increases the chance^ of 
mutual cancelling of rounding errors during arithmetical 
operations. Thus, on rounding to three decimals, 

Q-5772 -> 0*577, 3-14159 3*142, 

0-0625 -> 0*062, 0-1875 -*>0-188. 

This procedure normally ensures that the rounding error 
will not exceed half a unit in the last digit retained. 

Rounding errors can accumulate in arithmetical opera¬ 
tions, and it may be necessary to carry an extra decimal 
digit, otherwise unnecessary, rounding off at the end of a 
chain of operations. The probable error will be consider¬ 
ably less than the maximum error, and can be estimated 
statistically. For example, the maximum error in the sum 
of n rounded-off values is 'jzin units in the last digit and 
occurs only when each value has the greatest possible 
rounding error and ail reinforce; in practice the distribu¬ 
tion of rounding errors is likely to be almost random and 
for large n the probable error in the sum is likely to be less 
than ±4Vn/3, with error about ±Vit/3 in roughly one 
entry in twenty. 
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M 

Frequently we obtain finite numerical processes by suit 
ably truncating an appropriate infinite series. This neglect 
of higher terms introduces a truncation error , which can be 
kept small by retaining a sufficient number of terms 
(though with rapidly increasing labour as the number of 
terms increases). Truncation errors will be present in most 
processes other than those of arithmetic. 

Blunders are lapses from accurate computation due to 
carelessness of the computer or failure of the machine. In 
principle they are avoidable, but in practice they are almost 
impossible to eliminate; hence we must be constantly on 
guard, and the checking and recheckmg of calculations 
must become a part of our habitual practice as computers. 
For example, in adding a column of figures we should 
always do it twice, once from the bottom and once from the 
top; if the sums do not agree we must repeat until two suc¬ 
cessive checks do give the same answer. Whenever possible, 
checks should be made using different methods of calcula¬ 
tion, because if we have made a mistake once we are quite 
likely to make it again in a repeat calculation; careful com¬ 
putational design will often produce formulae which are 
not used directly in the calculation but must nevertheless 
be satisfied, and these provide valuable checks. We should 
always finish a problem with two consecutive checks that 
agree. 

Beginners often try to laugh off this emphasis on blun¬ 
ders, but they learn how easy it is: to transpose digits, 
taking 7239 in place of 7329; to misread repeated digits, 
taking 72239 in place of 72339; to misread tables; and to 
overlook signs, especially near sign changes. They find that 
after identifying a mistake they will make a wrong correc¬ 
tion. These and other blunders can be kept to a minimum 
only when the computer learns to work with an easy 
rhythm, to write his figures clearly and to lay out his com¬ 
putation in an orderly fashion. And he will discover his 
remaining blunders quickly if he checks the work stage by 
stage, avoiding checks by repetition; finally he will be con- 

11 
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fident of his answer only when he has completed an over¬ 
all check. 

Inexperienced computers often carry too many—some¬ 
times far too many—decimals through their calculations 
from lack of confidence, and from a vague idea that they 
may increase the accuracy of their work. We must under¬ 
stand clearly that although our model theories in science are 
capable of unlimited accuracy mathematically, they are 
based on approximate physical assumptions and are directed 
towards approximate physical measurements. In many ways 
the bridge of mathematics in between is expendable, and 
there is certainly no merit to be gained from unreasonable 
pursuit of accuracy in our model theories. Indeed the very 
striving for excessive accuracy defeats its own purpose; for 
the computer, burdened by excessive length in his calcula¬ 
tions, will quickly tire and err. We shall do well to recognize 
as the mark of the experienced computer the fact that he 
plans his problem so that at each stage he carries just that 
accuracy required locally to produce a final answer as 
accurate as the experimental measurements and no more so. 

1,5* Effective numerical computation 

The criterion by which we must judge numerical pro¬ 
cesses is effectiveness . This involves mathematical considera¬ 
tions such as accuracy and availability of independent 
checks, but also what may best be described as economic 
considerations such as the time required to obtain specified 
accuracy. When there is a long run of similar calculations 
wc shall spend ample time in reducing the problem to the 
best possible form for computation, so that each case may 
be handled quickly and accurately with an over-all saving 
in time and effort. But if we have a single case to compute 
we shall use any already-familiar process, and the calcula¬ 
tion will be finished in less time than a preparatory reduc¬ 
tion. We shall, of course, use mathematical tables, slide 
rules, hand machines, and above all else high-speed elec¬ 
tronic computers whenever they will be helpful* 
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A glance at some of the larger texts will show that there 
are very many processes available for each type of calcula¬ 
tion, such as interpolation or integration, for example. No 
doubt the beginner would wish to see one process of each 
type labelled best, but this cannot be done. There are no 
*best processes '; the most effective processes will always depend 
on the particular problem in hand . There is great scope for 
the use of judgment in numerical analysis in such matters 
as the selection of an appropriate process, the detailed 
design of a calculation, the choice of working interval, 
number of working digits, and so on. Indeed it is this de¬ 
mand on the capacity for skilled judgment by the computer 
that provides the challenge and much of the interest of 
numerical work. 

In this short book and a companion volume on Numerical 
Solution of Equations we shall discuss the kinds of way in 
which numerical processes can be constructed and put to 
use. We shall introduce few out of the vast array of avail¬ 
able formulae, but as the beginner enlarges his experience 
he will Learn how to select the most suitable from among 
many processes that might be used on a given problem. 


EXERCISES ON CHAPTER ONE 

1. State the number of decimal places and of significant figures in 
each of the following numbers, and (approximately) the accuracy to 
which they are presented: (i) 980*665, (ii) 911x10 (Iii) 2*9978 XlG 10 , 
(iv) 1-20x10-*. 


2. Compare the result of rounding off 1*37455 digit by digit until 
only two decimals remain, and rounding off to two decimals in a single 
operation. 

3. Evaluate and write In appropriate form, giving an estimate of the 
maximum error in each answer: 


(i) 1*73-2*16+0-08+1*00—2-23 -0-97 +3*02, 

(ii) 9*37x8-91, 

(iii) 7-137/0-279, 

(iv) 2*9132+15*164-0*31725. 


4. Compare the sum and greatest rounding error when the following 
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terms are added and then rounded to six decimals with those obtained 
when the terms are first rounded off and then added: 


*5833133 +*1000000 + *0083333 + 


*0011905+ *0002226 + 

+*0000493+*0000122+*0000033, 


5. Obtain linear, quadratic and cubic polynomial approximations for 
e* near *=0 by truncating the Taylor series appropriately, and use the 
remainder term to estimate (to the nearest 0T) the range over which each 
polynomial can be used if the error of approximation is not to exceed 
0*005 (he, the approximation is to be correct to two decimals)* 
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CHAPTER TWO 


Finite Differences 


Before we can develop particular processes we must in¬ 
vestigate the operations that can be carried out directly on 
a table of values. Continuous functions must normally be 
represented numerically by tables of values at specified 
intervals of the argument, and the limit processes of in¬ 
finitesimal analysis can no longer be applied to such tabu¬ 
lated functions but must be replaced by equivalent finite 
processes, 

2«i* Tables of values 

Most functions in common use are tabulated, and these 
tables of values provide one of our sources of raw material. 
Any scientist should become familiar with selected sets of 
tables: for example Comrie, Chambers's Shorter Six-figure 
Mathematical Tables is good for ordinary functions, and 
Jahnke & Emde, Tables of Functions is useful for some less 
common ones. To illustrate table notation we may take as 
typical Chambers’s table of ordinary logarithms 1,000 (1) 
10,000, 6D; the table gives in a rectangular array the 
mantissae correct to six decimals (denoted by 6D) for the 
numbers from 1,000 to 10,000 (the endpoints and the tabu¬ 
lar interval 1 in the argument are denoted by 1,000 (1) 
10 ,000); decimal points and characteristics are suppressed 
as they convey no additional information of value, and their 
inclusion would make the table harder to work from. Tabu¬ 
lated values correct to six significant figures are denoted 
by 6S* Another illustration is provided by Chambers’s 
table for e x given at values of the argument 0*000 (*001) 
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3 000 (■01) 10 00 (-1) 50 0 to between six and seven signifi¬ 
cant figures; the entries are in decimal form up to ar= 15*0, 
and thereafter in floating point decimal form 6D, 7S so 
that the entry for e ls =3,269,017 is printed as 3-269017, 6. 

Some functions change rapidly in isolated regions of the 
argument and slowly elsewhere, and such functions can 
best be tabulated by using different tabular intervals in 
different parts of the range according to the local behaviour 
of the function. The resulting tables are rather more diffi¬ 
cult to work from than those with uniform interval, and the 
great majority of the tables in common use are based on a 
constant tabular interval in the argument, at least over 
substantial regions (e.g. e x above). A section of such a table, 
say for /(*) = !/*, might be written 

* 30 3 1 3-2 3-3 3-4 3-5 3-6 3-7 

1 /* -33333 -32258 -31250 -30303 -29412 -28571 -27778 -27027 

This table cannot represent the nominal function 1/* 
exactly; for there will be rounding errors in most entries, 
and indeed the tabulated function is not formally defined 
at all for intermediate values of x. In this case we know very 
well that the tabulated values provide a good approxima¬ 
tion to ljx, and that intermediate values can be estimated 
by fitting a polynomial (Weierstrass’s theorem I) and we 
might anticipate (or hope?) that a low degree polynomial 
would be adequate because of the slow variation of l/x in 
this interval. But what if we are merely presented with a 
table of values? In principle this could represent some patho¬ 
logical function of pure mathematics, and even if Weier¬ 
strass’s theorem should be applicable a very high degree 
polynomial might be required for matching; or it might 
represent some perfectly well-behaved periodic function 
with natural wavelength considerably shorter than the 
range of the table. Thus, table 1 gives values for sin 2nx 
at three different working intervals h = 0-9, TO, and 1-05 
for the variable ac; it is ‘perfectly clear’ that each of these 
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Table 1 

* 0 A 2h 3h 4Ji bk 

f 0-9 0 —-588 —*951 —951 --588 0 

sin 2».« h= 1-0 0 0 0 0 0 0 

l 1*06 0 -309 -588 -809 -951 1-000 

three sets of values defines a ‘good smooth function', but 
unfortunately none of these ‘good smooth functions’ bears 
the slightest useful relationship to the original function 
sin 2irxl In each of these (rather artificial) cases the work¬ 
ing interval is so large relative to the periodic length of the 
function that the main periodicity has entirely escaped 
representation in the table. An even worse example is pro¬ 
vided by the function sin (2rrlx), for in this case as * de¬ 
creases, no matter how small a working interval we may 
select it will soon become impossibly large relative to the 
continuously shortening periodic length of the function. 
Clearly there can be no certain way of knowing whether a 
given table provides a proper representation of the (un¬ 
known) function which it is supposed to represent. The 
interval adopted in printed tables will always be short 
enough to display fully the essential character of the func¬ 
tion, and it is important that we should always adopt the 
same convention ourselves when giving tables of values; 
and it should be noted that this will usually mean using 
different intervals in different parts of the range according 
to the local behaviour of the function. We shall assume in 
working from tables that the higher order derivatives of 
the function are small enough for the function to behave 
‘reasonably smoothly’ over the length of the tabular 
interval. 

2.2. Finite differences 

, The simplest operation on a table of values, obtained 
either for a known function or by calculation, is to find the 
difference between each pair of tabular values. The first 
differences are found by subtracting each value from its 
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successor in the table, second differences by repeating a 
similar operation on the first differences, and so on for 
higher orders; these together comprise the finite differences 
of the table. The customary form for tabulation is: 


X 

/(*)= 

first- 

second- 

third- 


1/* 

differences 


3-0 

•33333 

-1075 



3-1 

•32258 

-1008 

67 

-6 

3-2 

*31250 

-947 

61 

-5 

3-3 

•30303 

-891 

56 

-6 

3-4 

•29412 

-841 

50 

— 2 

3*5 

■28571 

-793 

48 

-6 

3*6 

•27778 

-751 

42 


3*7 

*27027 





Note: (1) that odd-order differences are written between, 
and even-order differences on the tabular lines; (ii) that 
decimal points and leading aeros are omitted in the dif¬ 
ferences in order to simplify the table; (iii) that we must 
very definitely not omit the signs; (iv) that each difference 
must be multiplied by ICh 5 before use. 

Two features stand out in this table, (1) Differences of 
increasing order decrease quite rapidly in magnitude. We 
shall show in §2,4 that this type of behaviour is typical of 
polynomials and of functions or tables which are well-fitted 
for polynomial approximation. (2) The third differences 
show marked irregularity. If, however, we tabulate ljx to 
six decimals (6D) we find that the retention of the additional 
decimal has produced fairly regular third differences, but 
markedly irregular fourth differences: 
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2.2 


X 

1/* 

first- 

second- 

third- 

fourth 




differences 

3*0 

•333333 

-10752 




3*1 

■322581 

-10081 

671 

-60 


3*2 

*312500 

-9470 

611 

-53 

7 

3*3 

•303030 

-8912 

558 

-50 

3 

3*4 

■294118 

-8404 

508 

-40 

10 

3*5 

•285714 

-7936 

468 

-40 

0 

3*6 

•277778 

-7508 

428 



3*7 

■270270 






This behaviour can be seen to arise from rounding errors; 
even though individual rounding errors in /(jc) are ^i in 
the last decimal digit it is dear that these errors are likely 
to accumulate, and be increasingly amplified in succes¬ 
sively higher differences. The greatest error that can be 
obtained is, approximately. 


tabular 




error 

1st 

2nd 

3rd 

+i 

-1 



-* 

+ 1 

+ 2 

— 4 

+ i 

-1 

-2 

+ 4 

-i 

+1 

+ 2 

-4 

+* 

-1 

-2 

+ 4 

-i 

+ 1 

+ 2 

-4 

+1 

-1 

-2 



-i 


differences 

4th 5 th 6 tk 7 th 


+ 8 

-16 

-8 

+ 32 


+ 16 - 64 

+ 8 

-32 


-16 

-8 
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These are maximum errors, and the likely errors will be 
smaller as the rounding error will not usually alternate in 
sign; Comrie* has given a working criterion for the ex¬ 
pected fluctuations which will arise from round-off as 

order of difference 1 2 3 4 5 6 

expected limits of error ±1 ±2 ±3 ±6 ±12 ±22, 

and we must learn to tolerate small irregularities of this 
magnitude in our tables of differences. This level of fluctua¬ 
tions has been called the ‘noise level* of a difference table. 
The magnitude of the finite differences obviously de¬ 
pends on the length of the tabular interval. 

Although tables of values play a more important role in 
hand calculation than with high-speed machines, they are 
generally important for the expression of experimental 
results and of the results of a previous calculation. They are 
seldom used for the description of specified functions with 
fast machines, as it is usually much quicker to use a pro¬ 
gramme sub-routine to evaluate a function where it is 
needed than to feed in a table of numbers, 

2.3, Notation for finite differences 

There are three distinct notations for the single set of 
finite differences of a function f(x) tabulated at constant 
interval h in the independent variable x. If we write 
+ for the nth tabular point, and fn=f(x n ) t then 
we define 

(i) The forward difference at X{ v& fm—fu written Afu 
Thus 

^/(*o) + h) -f(x 0 ) =/l -/o, 

and higher differences A 2 fo = A(Af 0 ) — Afi — Aff =fz — 2fx +/q, 
j3/ £ = A {Affz) =/& - 3/4 + 3/3 —/2, etc,, can be calculated 
either by differencing lower differences or from the tabu¬ 
lated values. Note that 

• Sec introduction to Chamber*s Six-figure Mathematical Tables, 
VqL IL 
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2.3 NOTATION FOR FINITE DIFFERENCES 
M*}i =/* + , - ({)/*+,-> + (*) A +( - a - ... +(-!)*/*, 
where TM are the binomial coefficients. 


x 

x Q 

*1 

*2 

*3 



x n- 3 fn -3 

*n-Z fn-2 

x n -1 fn -1 

x n f n 


Vfn-2 

F%-i 



(ii) The backward difference at x% is ft—ft- 1, written 
pfu Thus, for example, 

V 2 fn = V{fn ~fn-l) ~fn ~ 2// t -l +fn-2< 

The forward and backward differences represent in 
a gg re § at e precisely the same set of numbers (and 
pft s but forward differences are specially useful 

near the start of a table as they involve only the tabular 
values below the ‘subscript level* Xi t and backward differ¬ 
ences are useful near the finish of a table where we do not 
have the forward tabular values to obtain all the necessary 
forward differences. 

(iii) The central difference at X{ is ft+%— written 
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S/|. Thus 8/(+Vs =/i+i ~fi = Afi = p/i+i, etc. This is the 
normal notation for use away from the ends of a table, 
where there are plenty of tabular values above and below 
the level *j. 


1 1 

>-*■ 14 

f 

5/(-% 


8/(+ l & 

/l+l 


*/*+»/. 

/1+2 


8 /ft- s /s 

fl+2 


SVm 

ayi 

8 2 /i+i 

8 a /i+a 


syi-vi 

S 3 /(+% 

S3 /<+v, 


sy« 

S 4 /m 


Note that d n /( is the «th forward difference and lies on 
the downwards diagonal from the tabular entry / (t while 
the nth backward difference p 7”f { lies on the upwards 
diagonal originating at ft. The wth order central difference 
S n / ( lies on the horizontal level x t \ thus odd central dif¬ 
ferences such as S/ifc lie on half tabular lines and even ones 
like S 2 /-i on tabular lines. 


2.4. Finite differences of a polynomial 

The finite differences of a general polynomial expression 
can be found only in terms of a specified tabular interval h. 

The finite difference of a constant over any interval is, of 
course, zero. 

The successive forward differences of the power x k are 

Ax* - (* + h) * - x*= khx *~i + *A(A - + ... 

+ AA* _1 *+A*, 

AH* - A(Ax*) - k(k - l)hH*~»+k(k - l)(A-2 

+ ... +{2*-2)k*, 

••*•!** 
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(k-t)¥^x*-^+ .... 


A*x*=k\ k*. 

All higher differences A* H x*, i>0 are zero. Thus the «th 
difference of x* is a polynomial of degree k—i for 0 <i ^ k, 
and is zero for t>k. 

The finite differences of an »th degree polynomial 

f(x)=ao+aiX+aiX 2 + ... + a n x n 

follow, and it can be seen that Af(x) is a polynomial of 
lower degree n-i if 0 <i ^ n, and is zero if t>n. In par¬ 
ticular, 

A n (fl 0 + &ix+ ... + a n x tt ) = a„n ! h n , 

d n+1 (<fO + tfl* + ... +«n* B )=0. 

Thus the nth difference of a polynomial of degree nisa constant 
proportional to h n , and higher differences are zero. 

Example 1. Derive finite difference tables far the polynomial 
/{x)=2—3x z +.'v 3 with the two working intervals A=0-6 and 

A = 0-1. 

These tables can easily be obtained as: 


X 

/<*) 

df 

A 2 f 

Ay 

Ay 

0 

2-000 

-626 




■B 

1-375 

-1375 

-750 

760 


1-0 

0 

-1375 

0 

760 

0 

1-5 

- 1-375 

— G 25 

+ 760 

750 

0 

8-0 

- 2-000 

+ 876 

+ 1600 

f 

V 

2-5 

- 1-125 


■' 


* 


* » 


* t * 
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X 

/(■*) 

Af 

A 

A*f 

A*f 

0 

2-000 

-29 




■1 

1-971 

-83 

-54 

6 


•2 

1-888 

-131 

-48 

6 

0 

•3 

1-767 

-173 

-42 

6 

0 

*4 

1-584 

-209 

-36 

* 

* 


-5 1-375 

• 

Note: (i) that in each case third differences are constant and 
equal to 3! and fourth differences are zero; (ii) that the dif¬ 
ferences depend on h and higher differences are systematically 
small only when h is small; (iii) that if tabular values of /(*) 
are rounded off there will be small irregularities in fourth and 
higher differences. 

A comparison of the second table of example 1 with the 
table for 1/x on page 18 suggests that after allowing for 
round-off fluctuation we might very well adopt a cubic 
approximation for 1/jc in the range 3*0 ^ 3-7. Indeed we 

shall find that whenever the higher differences of a table 
become increasingly small , the function represented by the 
table is well-fitted for polynomial approximation . Large dif¬ 
ferences may be reduced by reduction of the working in¬ 
terval hj but in some cases approximation by reasonably 
low-degree polynomials will be possible only for incon¬ 
veniently small h and other types of approximation may be 
more effective. Thus* it will be found that the values in 
table 1 are ill-suitcd for polynomial approximation and will 
remain so unless the interval is reduced considerably; this 
is hardly surprising for a sinusoidal function of wavelength 
small relative to the range of and we might expect that 
harmonic approximation would be better* 
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2.5* Detection of mistakes by differencing 

Rounding errors cause only small fluctuations in dif¬ 
ferences, but blunders will produce much larger upsets. 
Hence a standard and invaluable method for detecting 
mistakes in any table is to form differences and inspect 
those of higher order for irregularity. An error of one unit 
in f{x) (or in a difference)* and a pair of adjacent errors 


propagate through the differences 

as 

follows: 



0 




1 

0 




1 

0 


0 


1 


0 


1 


0 

0 


1 

-6 

0 


1 


-3 

0 


1 


-5 


1 


-2 


0 

1 


-4 

15 

1 


-1 


2 

1 


-3 


10 


0 


0 


1 

-2 


6 

-20 

1 


-1 


2 

-1 


3 


-10 


-1 


2 


0 

1 


-4 

15 

0 


1 


-3 

0 


-1 


5 


0 


-1 


0 

0 


1 

-6 

0 


0 


1 etc, 

0 


0 


-1 


0 




0 

0 


0 

1 etc. 

0 






Note how in each case the fanwise spread of the error 
locates its position, and that a single error grows on propa¬ 
gation as the binomial coefficients of increasing order. Also 
that the sum of the errors in each column crossing an error 
fan is zero. Hence, if we can establish the general trend of 
the affected differences in a column from those above and 
below an error fan* we can infer what these differences ought 
to be, Le. what correction is needed. In the simplest case 
of approximately constant differences the errors are ob¬ 
tained by a comparison of mean and actual values of the 
differences. 

Thus we detect errors in a table of values from the fact 
that a single error e produces maximum errors 

€ 2e 3c 6c 10c 20c 

in the 1st 2nd 3rd 4th 5th 6th dif- 
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DETECTION OF MISTAKES 


ferences, respectively. This is a powerful check and should 
always be applied to any table produced by calculation . 



of about 60 in 4 / 1*7 produces reasonably uniform fourth dif¬ 
ferences, and on checking we find that d/i*? should in fact be 
5716* 


Example 2* Check the following table for log^x. 
X log x differences 

1-0 '00000 

9531 

Irl '09531 -630 


recalculated differences 
9531 


-830 


8701 


133 8701 


133 



2'0 *69316 

The third differences indicate mistakes, probably more than 
one* Apparently and ^ 3 /i*3 are high, d 3 /i *2 and d 3 /i*4 low, 
suggesting too large an entry for log ff 1*4; the necessary reduc¬ 
tion lies between 20 and 30 (as may be verified by rcdifferenc- 
ing), and this suggests very strongly that the digits 74 have been 
interchanged from a correct value log e 1*4 — ■33647- The dif¬ 
ferences have been recalculated after this correction, and the 
surviving error fan is apparently centred on 4 / 1 . 7 ; a reduction 
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Note: (i) it should be verified that the fourth differences are 
now as uniform as can be expected; (ii) a convenient check can 
be made by comparing the sum of entries in any column of 
differences with the difference between the first and last entries 
of the preceding column—this check discloses the existence but 
not the position of the mistake ; (iii) differencing will uncover 
mistakes but not systematic errors* 

Problem 1. Tabulate s 3 —2* 2 -b 3*-1 to 4D for x=0(*Q 1)0*1, 

We could write a programme to calculate and print out the 


X 

/(*) 

F 

F 2 

F 3 

check 


-•02 

-1-060808 




- 1-0608 





30607 




306 


-■01 

-1-030201 


-406 


-1-0302 


-4 



30201 


6 


302 


•00 

-1 


-400 


-1 


-4 



29801 


6 


298 


•01 

-0-970199 


-394 

1 1 

-0-9702 


-4 



29407 1 


", 6 ’ 


294 


•02 

- -940792 1 


-38S* 

- -9408 


-4 



,29019*' 


6 


290 


■03 

- -911773 


-382 


- -9118 


-3 



28637 


6 


287 


■04 

- -883136 


-376 


- -8831 


-5 



28261 


6 


282 


•05 

- -854875 


-370 


- -8549 


-3 



27891 


6 


279 


•06 

- -826984 


-364 


- -8270 


-4 



27527 


6 


275 


■07 

- -799457 


-358 


- -7995 


-3 



27169 


6 


272 


•08 

- -772288 


-352 


- -7723 


-4 



26817 


6 


268 


•09 

- -745471 


-346 


- -7455 


-3 



26471 


6 


265 


•10 

- -719000 


-340 


- -7190 
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value of the polynomial at each successive value of x if we were 
using a high-speed machine, but this approach is unnecessarily 
laborious if we are working by hand. For we can build up suc¬ 
cessive values of the polynomial by accumulation of differences 
if we start from a few calculated values and use the fact that 
— const. 

The section of table above the stepped line has been obtained 
by summing the terms of the polynomial, using values sym¬ 
metrically disposed about x=Q for simpler calculation. The 
value/(*03) is then obtained from the sequence of operations: 
rV(-03) = -000006, rV(-03) = f 2 /(-02) + r 3 /<-03), f/(-03) 
= F/(‘02) + f 2 /(-03), /(-03) - /(-02) + p/(-03). This is in- 
dicated by arrows. 


Nate: (i) that even though the third difference is considerably 
smaller than the rounding error in four decimals it cannot be 
neglected because this would constitute a systematic error which 
would accumulate; (ii) that there is an easy check by substitu¬ 
tion at (iii) that the rounded values (which make up the 

final solution) have been checked by differencing; (iv) that 
backward differences are used as these depend on values 
already calculated. 


2*6* Divided differences 

Suppose that we have to deal with a function /(*) tabu¬ 
lated at values Xt of the argument x f where the tabular 
intervals are unequal . This calls for a more general type of 
finite difference taking account both of the change in tabu¬ 
lar value fi+i —/* and of the length of interval #( +l — xt; 
the obvious choice is the ratio (ft±i-ft)l(x t +i-xt) which 
gives a measure of the change in tabular value per unit in¬ 
terval. We define the first divided difference as 

[*,*♦,] . *”«■> -/(*■> . 

Xt+l-Xt 

Similarly the first divided difference for any pair of tabular 
values ft and f, (not necessarily consecutive or positively 
ordered) is 
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*/—*1 

The second divided difference [x ( x< + i.v (+2 ] of the tabular 
values f(xt), f(xt+i), f(xt+i) defined as 

r„ _ _ -i [*l+l#l+s]-[*<*l+l3 

[**«*«]-- letrt—xt -’ 

the third divided difference as 

Xi+Z-Xi 

and so on for higher orders. 

The divided difference of a set of adjacent tabular values 
occupies a position in the difference table at the apex of 
the triangle of diagonals based on the tabular values, and 
its value and position are entirely independent of the order¬ 
ing of terms in the divided difference. This may be seen if 
the divided differences are written symmetrically in terms 
of the tabular values ; for example, 

[x$xiX 2 ] 

_ [xix 2 ] - [xqXi] _ 1 f /(*a) -f(x 1 ) /( xi) -/fro) ) 

X 2 "Xq X 2 — Xq\ X 2 — X 1 X\—Xq f 

- /(**) , /frl) , /(**) 

fro - Xl)(xo - Xt) {xi - Xo){Xl - Xt) (xt - Xo)(X 2 - Xi) 

Similar results hold for differences of all orders (proof by 
induction). It can be seen from the symmetric form of these 
expressions that the ordering of the arguments is not 
significant, and that the value of the divided difference 
depends only upon which arguments are present; thus, for 
example, [xopcaxi;V 2 \ = [pcqxiX 2 X$]* 

The nth divided differences of a polynomial of degree n 
are constant and higher differences zero; for if f(x) = x* f 
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M = [arc-*] = — ' %Ll = *»-! +a;i»-a A - 0+ ... + . v j-i 

X — Xo 

so that the first divided difference is a polynomial of degree 
n — 1, and the result follows by repeated differencing. 

Example 3, Draw up a table of divided differences for 
tabulated at *=0, 5, 9, 12, 14, 15. 


X 

0 

/(*) 

0 

divided differences 

25 

6 

125 

151 

14 

1 

9 

729 

333 

26 

1 

12 

1728 

508 

35 

1 

14 

15 

2744 

3375 

631 

41 


0 

0 


The divided differences for a table of constant tabular 
interoath arc simply related to the finite differences intro¬ 
duced in §2.3; in particular, the nth divided difference of 
n+l consecutive tabular values is equal to A n fj(n\h n ) (and 
it should be noted that A*f l =p*f t+n = 8 n / (+B/2 , as these are 
alternative names for a single difference). The proof is by 
induction. Suppose that [wm x i+z ... x i+n ] = A n f i l(n\ k n ) 
for all i; then 


[*,*,+1 ...* (+n+1 ] - - [*«*<« ■•-*«+»] 
= {A n fMl(n ! h n ) - d”f(l(n ! h n )}j(n + 1)A 


The result is clearly true for first-order divided differences, 
and hence it follows that it is true generally. 


Problem 2- Evaluate (x— a&X*™ xi) (*-**„)[#* 0 * 1*2 #„] 

at x-x h i'=0,1, ft* 
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When two of the arguments in a divided difference are equal 
the formal definition reduces to 0/0 and further consideration 
is necessary. However, 

{x—Xi)[xXo Xi ... x n ] 

^(x-xt)[xx o ... Xi-lXi+l x n x t ] 

_ - _ x J*Q ... Xj-iXj+1 X n Xj}- [jfigQ ... Xj-lXj+l ... x n ] 

= [Aq ... *1—1**141 ... *#] [Aq ... *f—l*i*(4l .** *tf], 

and hence as *-►#(, (*— xt)[xx 0 x t ... 


EXERCISES ON CHAPTER TWO 


Little attempt has been made to provide straightforward exercises 
this chapter as such exercises can be constructed with ease by anyone 
possessing a set of tables. Readers are urged to work such exercises so 
that they may become familiar with the various processes and the ideas 
underlying them. 


I, Show that: (i) A [7/f = f/Afi = , 

(ii) Arfi = p n ft+n = 

(Note that in each case these are merely alternative ways of denoting 
the number occupying a particular position in the difference table.) 


2. Show that; (i) A(fm)^ft4gi+8Ui 
(til) J (Jtlgi)=(g t dfi -Jidgi )/C?(£ j+l ), (iv) 2(1 lfi)*= 


3. Show that: (i) J’*(l/*)= 


(-!)"»! 6” 
x(x+h) ... (x+tJi)’ 


(ii) 

where k is the length of the working tabular interval. 


4. Find tlie polynomial of lowest degree which fits the data: 

* 0 1 2 3 4 5 

y 75 150 204 241 265 280. 


5. 

it. 


rind which tabular value is wrong in the following table and correct 


* jW 

1-0 -340735 

M -500429 

1-2 -662671 

1-3 -825289 

1-4 -986092 


X 

/(*) 

1*5 

1*142289 

1-6 

1*293533 

1*7 

1435910 

1-8 

1 567997 

1*9 

1-687869 


X 

Ax) 

2*0 

1-793724 

24 

1*883902 

2*2 

1-956910 

2*3 

2-011438 

2*4 

2 046371 
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CHAPTER THREE 
Interpolation 


We assume here that tables of values provide good repre¬ 
sentations of well-behaved functions. Hence when we seek 
to interpolate, or estimate non-tabular values of a function, 
we take for granted that the function will behave smoothly 
between tabular points and that reasonable approximations 
can be obtained as polynomials of moderate degree. A 
polynomial of degree k will have constant Ath differences; 
hence we can judge the suitability of a table for polynomial 
representation by the behaviour of its higher differences, 
and when the kth order differences are nearly constant the 
function behaves in this range very much as a polynomial of 
degree k. 

3.1. Linear interpolation 

When the first differences of a table are nearly constant 
the function which it represents must behave very nearly 
linearly and can be approximated in each interval by a 
straight line joining adjacent tabular points; this is the basis 
of the familiar process of linear interpolation used almost 
universally in tables of logarithms and trigonometrical 
functions. Thus in the interval xo<x<xi, 

+ £=£{/[» 

Al — XQ 

When we are using tables of constant interval A it will be 
convenient to write x-xo- 0(xt - xq) = 6h , where 0 < 8 < 1, 
and then in the interval (#*, jff+j) 

f( Xi + eh) + e(f m -f t ) = (i - d)fi +*/, +1; 

or, using the various (entirely equivalent) expressions for 
the finite difference. 
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3.I LINEAR INTERPOLATION 

f(xi + &h) —fi + 6Afi = /1 + 6vfi+i = /< -f 0S/<+v*. 
Linear interpolation provides a simple process but one that 
can only be successful when the first differences vary very 
slowly, and this implies close tabulation to an extent that 
will be practicable only for very commonly used functions 
(e.g* "log* tables). In most tables we need interpolating poly¬ 
nomials of higher degree and we can no longer base an 
estimate on just two tabular values* 

Example 4* The following extract is from Knott's 4D log 
tables. 


X 

0 

1 

2 

3 

4 

5 

6 

7 

10 

0000 

0043 

0086 

0128 

0170 

0212 

0253 

0294 

11 

12 

13 

0414 

0792 

1130 

0453 

0828 

0492 

0864 

0531 

0S99 

0569 

0934 

0607 

0645 

0682 


proportional parts 


8 

9 

1 

2 

3 

4 

5 

6 

7 

8 

9 

0334 

0374 

4 

8 

12 

17 

21 

25 

29 

33 

37 

0719 

0755 

4 

S 

11 

15 

19 

23 

26 

30 

34 


Differencing the values of the first row: 


X 

A*) 

A 

A* 

100 

0000 

43 


101 

0043 

43 

0 

102 

0086 

42 

-1 

103 

0128 

42 

0 

104 

0170 

42 

0 

105 

0212 

41 

-1 

106 

0253 

41 

0 

107 

0294 

40 

-1 

108 

0334 
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INTERPOLATION 


The differences are locally almost constant and the table is 
suitable for linear interpolation: thus 

/(1033) = 0128 + 3 x 42/10=0141. 

The subsidiary tabic of proportional parts shows tenths of the 
mean difference between tabular values for a row of entries; it 
should be noted that the values printed are a compromise, and 
that the proportional parts between 100 and 102 should be 
4 9 13 17 22 26 30 34 39 with those between 107 and 110 as 
4 8 12 16 20 24 28 32 36, and a single table can serve the 
range (100, 110) only at the cost of some loss in accuracy (very 
nearly of the fourth decimal digit!). Good tables do not use 
averaged proportional parts because of the consequent loss in 
accuracy. 

In contrast, if we difference the first column we obtain the 
table: 


* 

/(#) 

A 


A* 

10 

0000 

414 



11 

0414 

378 

-36 

5 

12 

0792 

347 

-31 

6 

13 

1139 

322 

-25 


14 

1461 





In this case we should require a polynomial of degree about 
three for interpolation, judging from the higher differences; 
hence this tabular interval is inconveniently large for a working 
4D log table. 

3.2. Interpolation using divided differences: New¬ 
ton’s formula 

Suppose that we have a table with constant tabular inter¬ 
val and we want to interpolate a value. The very act of 
putting in the ‘new tabular point’ at which we wish to 
determine/ destroys the original uniformity of the tabular 
interval by producing two smaller intervals about the new 
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3.2 newton’s formula 

point, of length Bh and (1 — 6)k, so that we must for the 
present work in terms of divided differences. 

The problem remains essentially the same if the tabular 
intervals are generally unequal as in §2.6, and we shall 
assume this for the time being. We now seek to represent 
the value /(*) at a non-tabular value x in terms of finite 
differences for the tabular points, and this can be done in 
the following way. The original definitions for divided dif¬ 
ferences of successively higher order can be rewritten as 

f(x) = /(# 0 )+ (# - x 0 )[xx 0 ], 

[**o] = [*0#l] + (* - #l)[##tt*l], 

[##0#l] = [#0*1*2]+ (* — *2)[**0*l*g]t 


... #n-l] = [*0*1 ...#«]+(# — *„)[**o*i 

By successive substitution, for [**o] in the first equation 
from the second, for [##o#i] in the second from the third, 
and so on, we obtain 

/(*) =/(*o) + (* - *o)[*o*i] + (# - * 0 )(* - #1) [*0*1*2] + .. . 

+ (*-#o)(#“#l) ... (#-*»-l)[#0*1*2 ...*»] + R n (x) 
= Pn(x) + R„(x), 
where 

Rn ( x ) = (* - #o)(# - #i)(*-*2)... (*-*„) [**0*1 ...*„]. 

It should be noted that *o, *t, ... , * n can be any selected 
ordering of the n given tabular points, and that * can lie 
anywhere in the general range of the tabular points (and 
in particular need not be greater than *o). 

This expression for /(*) would immediately provide an 
interpolation formula of the kind we are seeking (/{*) given 
in terms of tabular values and their divided differences) if 
only we could assert that the remainder term £„(*) was 
small relative to /(*). When the tabulated function is a 
&th degree polynomial, divided differences beyond the Ath 
are all zero and /?„(*) = 0 provided that n>k: hence we can 
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INTERPOLATION 


interpolate polynomials exactly (or, to be more precise, to 
within rounding error). 

Example 5. Find (9-S) 3 from the table of example 3. 

We can identify xq, *i, ... with the given tabular points in 
any order we please, and it will be sensible to select them in 
increasing order of distance from the non-tabular value 9*3 
so that the factors ( x —ay) in -P fl (j;) may be kept as small as 
possible. Hence take * 0 =9, *i=12, x z =5, x 3 = 14, *4=15, and 
# 5 = 0 ; and then 

/(9 P 3)= 729+ 0-3x 333+0*3x (- 2*7) x 26+ 0-3 x (-2*7) 
x 4*3 x 1 + 0-3 x (- 2*7) x 4*3 x (- 4*7) x 0 + 0 
= 804*357 

This result is, of course, exact (and also trivial, though it serves 
very well as a demonstration)* Normally we should work to a 
specified number of decimals and there would be rounding 
error in the table and in the final result* 

The functions we have to interpolate will not usually be 
polynomials, but whenever they are well suited to poly¬ 
nomial approximation the higher divided differences will 
be small (and indeed this is how we judge suitability). In 
such cases we use the interpolation expression without 
remainder, 

f( X ) ~f( X o) + (x- tfo)[#0#l] +(#“ Xo){x - Xi)[X qXiXz] + 

* "b #o)(# #l) »,*(# — Xft -l)[#0#1.• * • 

This is Newton's interpolation formula * The choice of order 
(n) for the approximation is a matter for the judgment of 
the computer and depends on the behaviour of the higher 
differences in the table* 

The truncated Newton interpolation formula 

f(x)^P n (x) 

is a polynomial approximation of degree n in which the 
coefficients are given in terms of the divided differences of 
the table. The remainder may be regarded as an 

error term; and provided that f(x) has n+1 derivatives in 
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newton's formula 


3‘ 2 

a range including xq , *i, ... , it can also be expressed as 

R n (x) = (x - xq)(x *»)/<*+«( £)/(n + 1)1, 

for at least one value i in the range (see problem 3 below). 
It may be seen from this form for the remainder or error 
term that for a given function we shall have a small in¬ 
terpolation error if the tabular intervals are all short and 
if we refrain from introducing tabular values distant from 
the level x. 

A considerable practical advantage of Newton*s formula 
is that we can at any stage add or utilize an additional value 
anywhere in the range without disturbing the terms already 
calculated. In general f(x ) ^P n (x) except at the tabular 
points, but the error R n (x) lies within assignable limits; and 
for small R n the tabular points x 0f ..., x n should be chosen 
as nearly as possible symmetric about x . P,*(#) is the 
smoothest function that agrees with f(x) at the tabular 
points, as d n+l P n (x)jdx n+l = 0 for all x* 

It may be worth emphasizing that the symbol ^ is used 
in the expression f(x) ^ P rt (x) to show that the right-hand 
side is a polynomial representation and that some question 
of judgment is involved in the selection of P n (x). We are 
asserting that the function is sufficiently well defined by the 
table of values and is well-represented by a polynomial of 
appropriate degree* Many authors use an ordinary equality 
sign, but we wish to draw attention to these steps of approxi¬ 
mation. 


Problem 3. Show that when / has n derivatives there is a 
relation of the form [xqx j# 8.„ #*,]=/<**>( |)/n I. 

Since j[x)=P n (x)+RJx) it follows that R n has n derivatives 
too. Hence, using Rollers theorem, as R n (x) vanishes at the 
? +1 points Xq , Xi $ ... Xn it follows that PJ,{x) has at least n zeros 
in the range, P£(#) at least n—l zeros, ,,, and ultimately that 
P^ J (#) has at least one zero, at some intermediate point £ of 
the range* But 

/(^W)+^ fl | [xoxt ... x n ]+Rp(x); 
and hence ,,, Xn]=f in) (£)[nl 
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INTERPOLATION 


3.3* Lagrange interpolation formula 

For some applications it is more convenient to have a 
matching interpolation polynomial with coefficients which 
are given directly in terms of the tabulated values of the 
function, as in the case of the Lagrange interpolation formula , 


f{x) ^ L n (x) 

„ {*-*iX*-**)»■ (*-*») f{x V . 
(Xq-Xi)(x 0 -X2) ... (xn-x n ) j{ 

(x~xfi(x-Xi) mm .(x-Xn) « , , 

(#—(#— Xn—x) ^ \ 
(Xn-Xo)(Xn-Xi) ... (X n -X n - l) A 


The polynomials L n and P n are identically equal. For 
L n is a polynomial of degree n which takes the values/(# ( ) at 
the tabular points Xt (* = 0,1, n ). Moreover, L n (x) - P n (x) 
is a polynomial of degree n with the n -f1 zeros Xt , and hence 
is identically zero;* thus L n is merely another form of the 
Newton polynomial 

In practice the Lagrange polynomial L n (x) is less con¬ 
venient for interpolation than P n (x)* The value f(x) must 
in general be more closely related to neighbouring tabular 
values ft than to distant ones; this special dependence on 
neighbouring values is recognized implicitly in the form 
of P n (x) and it is relatively easy to adjust the degree of the 
approximation to suit local conditions, white in L n (x) all 
arguments appear symmetrically so that all terms must be 
calculated and adjustments are less easy. 


3.4. Interpolation in a table with equal intervals 

Most interpolation will be in tables with constant tabular 
interval and we shall now develop from Newton’s formula 


* Recall that: just one polynomial of degree less than or equal to n can 
be fitted to » + l distinct points; a polynomial of degree n which vanishes 
it n-fl distinct points is identically zero* 
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3.4 INTERPOLATION WITH EQUAL INTERVALS 


other forms appropriate for use with constant interval 
tables. If we want the approximate value f(x) it is natural 
to choose Xq as the algebraically smaller adjacent tabular 
point to jc, so that x—xq +0h, where G<0<1 and h is the 
tabular interval; beyond this the choice of arguments need 
not be the natural order, and indeed we can obtain a variety 
of interpolation formulae by taking different orderings of 
the tabular points. 


Xq 

XQ 

h 

4/0 

4/i 


Xl 

XQ±k 

h 


*2 

xo+2k 

h 





4 h 


*3 

*0+3 h 

h 

4fr 


*4 

#o+ 4/i 

U 





A% 


The simplest case is that of interpolation near the head 
of a table, in which case we identify the points #q, Xi t ..* of 
Newton’s formula with the tabular points in increasing 
sequence (i.e. at progressively greater distance from #). We 
must use forward differences as these provide the only 
higher-order differences available near the head of a table, 
and hence Xt — Xo-tkh, x — xt = (6-k)h. Bearing in mind 
that [xqXi x n } — A n fol(u} h n ), we have by direct sub¬ 
stitution in Newton’s formula 


/(*)=/o+fj4/ r o + 


>-l) 


A*fa+- 




^ 3 /o + 




This is the Nezvton-Gregory forward difference interpolation 
formula for use near the start of a table, mainly in the 
interval (*o, Xi)> 

In a similar way we can construct a Newton-Gregory 
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INTERPOLATION 

backward difference interpolation formula; however, we 
shall be more interested in a Newton-Gregory backward 
difference extrapolation formula which may be used for 
extending a table beyond its last value* f(x n ). Thus if 
x=x n +Bh, we obtain 


/(*) —/n + ^ r/n + + ... + 


f?{l + 0) ...(«-1 + 0) 

+ ^- -7*U 


Most interpolation will be in the interior of tables, where 
we use central differences so as to give equal weight to the 
tabular values on either side of the working level. In this 
case we can introduce tabular values roughly symmetrically 
about the interpolation level if we identify the v 0 > *1 of 
Newton s formula with the tabular points in order' of 
increasing distance from x starting with the point alee- 
braicallv mat -_ _ F i 



xo~3k 

/- 8 



X4 

x 0 - 2 h 

f-2 

U 

S/-3/2 


X 2 

Vo— h 

8/-1/2 

sy_i 

fit 

*0 

Vo 

S/l/2 

sy 0 


*1 

vo+A 

v 0 +2A 

fl 

8 /3/2 

8% 

*3 

h 

sy 2 

vq+3A 

S /6/8 

*5 

h 


* 

vo+4A 

U 




/- 1/2 


S 3 /i/s 

S 3 /3/2 


5 4 /o 

8 4 /i 


illustrated. In this 


case x=xo+$k, x-x 2]c = (6+k)h, 

* Extrapolation is a potentially more dangerous pursuit than internola. 

Ty do d n« CaSl,y eXtrap0,ate funCti0n * into 
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3.4 INTERPOLATION WITH EQUAL INTERVALS 

*-*s*-i = {Q-k)h, and [vovi... x 2n ] = S 2 ”/o/(2«)! h 2n , 
[vovi ... * 2 t»+i] = 8 an+1 /t/)/(2n +1)! A 2rt+1 on account of the 
selection order of the tabular points; and on substituting 
in Newton’s formula 

f( x ) —fo+(x— * 0 ) [vovi] + (*—Vo)(v—* i)[j6q*ijc&] + ... 

we obtain the ‘basic Newton formula' for constant interval 
interpolation 

/(*) =/o+9S/*+fcl> 8y*+ 

+tL K^i) Msft+ _ 

Note that if we sometimes write our formulae as unter- 
minated series this does not mean that they are to he 
regarded in any sense as infinite series, but only that the 
truncation point is left to the discretion of the computer in 
each case; if polynomial representation is suitable, the error 
will at first decrease as we take extra terms, but beyond a 
certain stage it will generally increase again. The basic 
Newton formula can be rearranged in more useful forms. 
(1) For interpolation over the inner part of the interval 
a - o.Vi we can derive a form symmetrical about the mid¬ 
point a-o+J/i by rewriting even differences as, e.g. 

sy 0 =i(sy 0 +syo+*{sy 0 - syi)=i(sy 0 +syo - is®/*. 

Thus 

J (*) — {i(/o +/i) + if/o —/1)} + Bhfh + 

{i(syo+syi)+Ksya-sy 1 )}+... 

-K/c +/i)+(0-»s/. A +(sy 0 +syo+ 

+ ( sy 0 +syo + ..., 
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INTERPOLATION 


This is Bessel’s interpolation formula. Bessel’s interpolation 
formula is sometimes written 

/(*) — pf’h + (9 -$)Sfih + —Yi 
+ ... 

where the averaging symbol p denotes the mean of two 
adjacent elements of a vertical column; thus 2 pf‘h=fo+fu 
etc. It is also written 

/(*) =/o+ es/vt+ £*(sy 0 +s 2 / : )+ sr&fa + B iv (h%+ syi) 
+ 

where BT= ±6(6-1)121, BT = 6(6- *)(0-l)/3I, ... are the 
Bessel coefficients. 

(2) Some sets of tables include auxiliary sub-tables of even- 
order central differences (as these lie on tabular lines); to 
obtain an appropriate interpolation formula we need only 
replace the odd-order differences in the basic Newton 
formula by the one-Iower-order even-difference combin¬ 
ations, 


/(*) —fa + *(/i -fa) + m 

“I* • a * * 

Taking j> = 1 — 0 , 

, /o + J 

+*(a+V*+ S2 = 1 F=®**+™}. 

This is Everett’s interpolation formula, 

(3) We can obtain a form symmetrical about tabular 
points as follows: 
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3.4 INTERPOLATION WITH EQUAL INTERVALS 

fax) —yb+ o(8fk - isy„)+~ 1 w-fa i)(s»/* - isy 0 ) 

+1 Wo+... 

= fa + i^(S/+ S/'/s) 4- -1 ^S 2 /o + 

+~0(e 3 - 1 )(sy_.A+53/ 1/1 ) + ... 

=/o + OfiZf 5 + 2] ^ 2 fa + 3I ~ ^)f® 3 /o + 

-i)sy 0 +..., 

and this is tnter/>o/ai* 0 « formula \ we shall use this 

form later in deriving differentiation and integration pro¬ 
cesses. 

A useful formula for half-way interpolation follows from 
the Bessel or Everett forms by taking B — \: 

A ^ l[(fo +A) - + 3%)+j|g(sy 0 +sy,) - 

- j^w*+m}+o(h*). 

This is simple to apply and useful as a first step in breaking 
down a table with large interval. Note that this form of 
expression indicates approximation and also shows (in 
terms of the interval) the order of magnitude of the error 
term. 


Example 6. Use BesscFs interpolation formula to find 
illO/3) from the table for /(#)= 1/x given on page 19- 

Here 0=(lO/3— 3 * 3 )/( 3 " 4 — 3 * 3 )= 1 / 3 , so that the coefficients o f 
the Bessel formula are ${&— 1)/2.2!*= —1/18* 0 ( 0 —1)/31 
-1/162, {0+1)0(0- 1X0- 2)/2.4!=5/486,... Thus/ 
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INTERPOLATION 


/(3°) ^ *303030+1(- -008912)- ^{-000558+ -000508)+ 

+ li^<—’000050)+ ... 

=-300000^’ 00297 ° 7 ~ ' 0000592 ~ -0000003+ ... 
i.e./(3J) -300000 on rounding off. 

is not. J J t worth retai nmg, that from S 4 / 

3 ‘ 5 ‘ t£owb S acb° f BCSSel a ” d Everett * nter P°iation: 

.J[* v >! S0lated inter P°Iation we shall use the formula 
with which we are most familiar, but for a run of inter 

formX 8 " WU be advanta g e °us to select the ‘best’ 
Bessel's formula is one of the most generally useful T„ 

™rsirf fc, ulficient ,o work 

Oiuy (.giving a first correction to linear intemolarinn’i 
have 8 t s ° met!mes t}lird and occasionally higher differences 

of ST r, account - Wc ^ a wXI 
WKS35el5Sif 33 BKsd & 

d£ ^ d SS 

can be WH f t ° f the crror in an y Particular case 

but these rou ^ h 

Mathematical tables often include some finite differences 
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3,5 BESSEL AND EVERETT INTERPOLATION 

If only first differences are given the table will take linear 
interpolation; but when higher differences are needed it is 
common to print only 8 z f and S 4 / because of the difficulty 
involved in printing odd differences on half-lines, and in 
this case Everett’s formula will obviously be more con¬ 
venient. A further advantage is that second difference 
terms in Everett’s formula take complete account of S 3 /, 
and similarly Everett up to fourth differences actually in¬ 
corporates the effects of the first five differences; hence it 
will be superior when third differences cannot be neglected 
in Bessel’s formula (8 3 />60) but fourth differences are 
negligible, or when fourth and higher differences are sig¬ 
nificant (as when many decimal figures are needed, or 
when the tabular interval is large). It can easily be shown 
that Sy ^ 4 and 8 s / ^ 20 can be neglected in applying 
Everett’s formula with a consequent error of less than half 
a unit. 

In an effective interpolation formula we hope that there 
will be few terms to calculate and yet we must avoid exces¬ 
sive truncation error, and we can make gains in both these 
directions by taking partial account of one of the higher 
differences in tire procedure known as throwback. For 
illustration, in Bessel’s formula the whole contribution 
from the second and fourth differences is 


'■”fer > ( (s v»+ s %) - jg(i+e(2-w+«*/o). 


By a fortunate numerical chance the coefficient 
(l + 0)(2 —0)/12 varies little over the range 0<0<1; in¬ 
deed it lies between -1800 and -1875 for -2<0<-8, and 
near the ends of the 3-range the whole contribution is 
small because of the outer factor 8(8— l)/4. Hence we get 
a very good approximation to a Bessel formula carried to 
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Fig. 1. 


3,5 BESSEL AND EVERETT INTERPOLATION 

/(*)-( 1 - $)f Q + 0f :l - mO-1)(&~ + 

+$(e+me-i)hif u 

Summarizing: (i) near the start or finish of a table use a 
Gregory formula; (ii) when 8 3 /<60 units in the last digit 
use Bessel up to second differences ; (Hi) when S 3 />6G but 
S 4 /< 1,000 units (or 400 units if S 4 /o and S 4 fi are of opposite 
sign) use Everett with throwback from fourth to second 
differences; (iv) when S 4 /> 1,000 units use throwback with 
higher differences. 


Example 7. Find /(2/3) from the table for /= ifx given 
below* 

In the following table &=*05 and jc= 2/3= ‘65+h/S, so that 
6— 1/3* Third differences cannot be neglected, but as fourth 
differences are less than a thousand we should normally use 
Everett with throwback, and we need then take no further 
account of 8 4 /. As an illustration this particular example will be 
worked using: (i) Bessel to S 4 /; (ii) Bessel to 8 3 /with throwback; 
(iii) Everett to 8 %f. 


m 


S 2 


S3 & 


-C54 Sg, 


fourth differences if we work to third differences in the 
formula and use the modified second differences 

S2/< = S3/,-CS*/ ( , 

where C is a constant, and may be taken as 0-184; this is 
known as throwing back the fourth difference on the 
second. Some accuracy is lost (though this throwback error 
is unit for S 4 /= 1,000 units!), but we make a sub¬ 
stantial saving in labour, and this is precisely the sort of net 
balance we always seek in designing numerical processes. 

In the same way we can throw back fourth differences on 
the second in Everett’s formula, and again we use 8 m 2 /i 
e=S 2 /i —0-184 S 4 / ( in place of S 2 /j and the modified Everett 
formula is 
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*50 2 00000 

—18182 

*55 1-81818 3031 

-15151 -701 

•60 1-66667 2330 203 

-12821 -498 

■65 1-53846 1832 131 

-10989 -367 

'70 1-42857 1465 93 

- 9524 -274 

•75 1-33333 1191 63 

- 8333 -211 

•80 1-25000 980 

~ 7353 

•85 1-17647 


-37 2293 
-24 1808 

3256 

-17 1448 
-12 1179 
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Hence: (i) 

0((?-l)/4= —1/18 

1/162 

(0+ 1)0(0- 1 )(«- 2)/48= 5/486 

/o= 153846 

0S/ % « - 03663® 
-(8%+«yi)/18- —*001832 
S 3 /tJl62= — ‘00002a 
6{Sy 0 +Syi)/486= ‘000023 


{ii) 


/o= 1*63846 

8/W3=-03663o 

(BJW+MAVW- -‘001809 
8 3 /iyi62— -*000023 


1-499393 


l*49999g 

(iii) 

1 )/ 6 — -1/81 
fl(^-l)/ 6 = — 4/81 

/o= 1*53846 
3/J3= -*G3663o 
- 8 * 1 / 0 / 81 = - 001 lie 
-48j/i/81= —‘00071s 


1-49999® 


In each case the rounded solution is 1*50000- 

Note: (i) that the subscript m denotes a difference modified 
by throwback; (ii) that the effect of neglecting S 4 /is about 2 in 
the last place, and some account must be taken of S 4 / if the 
interpolate is to be as accurate as the tabular entries; (iii) that 
successive terms do not necessarily decrease mo no to ni call y in 
magnitude; (iv) that the probable round-off fluctuation of 16 
in S 4 /has negligible effect* 

3*6. Some uses of interpolation 

(i) Subtabulation. The choice of tabular interval in the 
construction of a table is influenced by the uses anticipated 
for the table, but also by the economics of book production; 
hence we shall sometimes find that the printed tabular 
values are too widely spaced in regions of special interest. 
The process of ‘filling in’ between tabular values (known 
in this context as pivotal values) by systematic interpolation 
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3*6 SOME USES OF INTERPOLATION 

is called subtabulation. For example, if the initial tabular 
interval h is to be broken down to A/ 10 , it may be worth 
while applying the halfway interpolation formula followed 
by Everett using throwback as necessary. Subtabulation 
should be carried at least to an interval such that accurate 
interpolation is obtained using no higher than second dif¬ 
ferences; once a table can be linearly interpolated there is 
no gain from further subtabulation* 

Note: (i) that if a table of a complicated function is to be 
built up by calculation from a formula and a fast computing 
machine is not available, the easiest procedure may be to cal¬ 
culate at large interval and subtabulate; (ii) that after halfway 
interpolation the differences in a table will be quite different 
and higher differences will be considerably reduced (as 2~ n ) t 
and that the usefulness of subtabulation arises from the fact 
that 8 ft /== 0 (A n ); (iii) that a table of irregularly spaced values 
can first be reduced to regular interval using Newton’s divided 
difference formula, and then subtabulated in the normal way. 

(ii) Inverse interpolation. In the process of interpolation 
wc seek the value of f(x) for a specified non-tabular value a?; 
the inverse process in which wc seek the value of x cor¬ 
responding to a specified non-tabular value f(x) is also im¬ 
portant (e.g. in extracting roots of the equation f(x) = 0 ), 
and this is called inverse interpolation. It can be carried out 
by successively closer subtabulation in progressively 
smaller neighbourhoods of the specified value f(x) until 
linear interpolation is possible (note that this approach 
allows full scope to the computer for an inspired guess at 
the neighbourhood and interval for starting). 

Inverse interpolation can also be carried out by methods 
of successive approximation; this approach cart be illus¬ 
trated using Everett’s formula. Everett’s formula to second 
differences can be rearranged in the form 

6 ~J^Jo[ f{x) ~ fo+ ~ ie+ Wd + ■••}• 

Provided that the function f(x) is well defined by its table 
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INTERPOLATION 

of values, the terms involving second and higher differences 
will normally decrease fairly rapidly, so that we can expect 
to get a modest first approximation 0i by ignoring the con* 
tributions of second and higher differences, and taking 

„ /(*)-/» /(*)-/« 

~ * 

To improve this approximation we must bring in the term 
in second differences, and this involves but as this term 
will provide only a relatively small correction to 0*, we can 
get a satisfactory second approximation 82 from 

0, = {/(*) -/ 0 + iB 1 (9 l -1 - 2)S*/o - (8 1 + l)Syi]}/S/% 

where the error introduced by writing 81 in place of 8 in the 
second-difference term is likely to be of the same order as 
the higher-difference terms which arc neglected. The 
accuracy of the approximation can be improved by the 
successive introduction of terms in higher-order differences 
(though only to the degree permitted by the quality of the 
table). Note that the number of decimals which should be 
carried at each stage can be estimated by making the round¬ 
ing error roughly comparable with the first term neglected 
(and it should be stressed that there is no gain from carry¬ 
ing too many figures in the early stages of the calculation)* 

Example 8. Find the value of x for which log e #*=0*3 from 


the following table. 






X 

l°gr r 

8 

8® 

8 s 

8< 

8® 

85 

1-2 

*182322 

80042 

-6970 

1036 

-292 

75 

-57 

1-3 

•262364 

74108 

-5934 

819 

-217 

59 

-16 

1-4 

■336472 

68993 

-5115 

661 

-158 

36 

-23 

1*5 

•405465 

64539 

-4454 

539 

-122 

33 

-3 

1-6 

■470004 

60624 

-3915 

450 

-89 

12 

-21 

1-7 

■530628 


-3465 


-77 


+ 10 
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3*6 SOME USES OF INTERPOLATION 


Working from the table, 

e = 0-3—/o _ -3--262364 

1 / 1-/0 -336472--262364 


•50785s. 


Our best choice is probably to proceed with the value 0* = Q’5, 
if we take due account of ease of calculation at the next step and 
also the likely error in this term! Then, using Everett up to 
modified second differences (we may just as well introduce at 
this stage the small but not negligible contribution from the 
fourth difference) we have for the second approximation. 


02= {/(*)-/o+ Wi- l)Pi- 2)3 %f 0 - ( e 1 + Wi]}/S/* 


= {•3-262364—![-§(- 


•49859s. 


-005894)—1( —*005086) 
-+*074108 


b 


A recalculation of $2 using the improved approximation 
0=^*49859 in place of 0^*5 will now serve both to reassure us 
that no substantial error has arisen from use of the crude 
value 0i=*5 in calculating 02 (he, that the approximation is in 
fact self-consistent to the second stage), and to check 0a against 
blunders* 

02={-3—262364+ Jx -49859 x (—50141)1-1-50141 x 

x (- -005894) -1-49859 x ( - -005086)]}/-074108 

= •498593- 

There will be no gain in accuracy by proceeding further, as the 
fourth differences (of about 200) have been sufficiently taken 
into account by throwback and the contribution to the numera¬ 
tor from the sixth differences is negligible. Hence 

* ±hl*3+ ‘49859 x 0-1= 1-349859. 

Note; (i) that the irregular appearance of the sixth differences 
is sufficiently explained by the probable round-off fluctuation 
of ^22; (ii) that by recalculating 0g we avoid carrying on an 
unnecessary error to the next stage, and in this case if we had 
used 0i= *508 we would have found 02—*49859 and might have 
carried on *49860; (iii) that processes of successive approxima¬ 
tion in which each improved estimate is based on the previous 
one are called iterative processes; (iv) that modified difference* 


















INTERPOLATION 

arc useful in inverse as in direct interpolation; (v) that the 
accuracy of the process is good provided that is not small, 
but that for small S/i^ the rounding errors can be disastrously 
amplified, for example near a stationary point of /(#), and it is 
then usually best to subtahulate until linear interpolation is 
possible; (vi) that an initial survey of the problem is absolutely 
essential so that the computer may have some idea of the sort 
of pitfalls that must be avoided; (vii) that when higher differences 
arc large it may be best to subtabulate at interval JA as a first 
step; (viii) that on no account should the interpolation formula 
be treated as a quadratic (or cubic) to be solved for 0, as this is 
a laborious and inaccurate approach. 


3*7. General notes on interpolation 

(i) Tables of Bessel and Everett coefficients for 0=0(*GGI)1 
can be found in Interpolation and Allied Tables (H,M*S.O.), 
in Chambers's Six Figure Mathematical Tables y VoL 2, and 
elsewhere. 

Critical tables of coefficients are specially useful, and a 
set of critical tables for Bessel coefficients is given in 
Chambers's Shorter Six Figure Mathematical Tables , Bessel 
coefficients vary slowly with 0, and it is possible in quite a 
short table to write down all the values taken by a co¬ 
efficient, for example 2T=‘001, '002, '003, .and to 
list beside them the critical values defining sub-intervals 
of 6 vriiich correspond to each value of B'\ The leading 
terms in the critical table for IT are: 


6 

*000 

*002 

•006 

■010 

■014 


B' 


— *000 
— *001 
— ■002 
— *003 


These critical tables are easy to work from and are very 
convenient, 

(ii) If instead of a table of values we start from a function 

52 


EXERCISES 


which can be differentiated with reasonable ease it may be 
more effective to evaluate the function and its first few 
derivatives at pivotal points of the range, and then to sub¬ 
tabulate using Taylor series based in turn on each successive 
pivotal point. The terms in the Taylor scries are in general 
smaller than corresponding terms in any finite difference 
interpolation formula, 

(ill) Scientific literature contains many accounts of (so- 
called) semi-empirical theories^ in which dimensional and 
other incomplete arguments are used to derive solution 
formulae containing several disposable constants, and these 
constants are to be evaluated from experimental results to 
complete the solution. But any given segment of curve can 
be matched to specified accuracy by a polynomial, which is 
fully determined by its coefficients. Put in another way, any 
curve segment (or set of experimental points) can be matched 
by a function containing a sufficient number of disposable 
constants, and the relatively simple curves which ought to 
be the product of wcll-conccived experiments can normally 
be matched over moderate ranges (and within experimental 
accuracy) by quite low degree polynomials. Hence we must 
treat semi-empirical theories with great reserve; for the 
link between theory and experiment has been reduced (at 
worst) to one of curve-matching and we have forfeited the 
direct relationship which alone justifies our construction of 
theories. 


EXERCISES ON CHAPTER THREE 


Again exercises may easily be constructed from any set of tables. 


1, Subtabulate the following 

table to produce 

a table of f{x) for 

* = 1*5(1/30)2-4* 


* 

1-5 

m 

*40546 

X 

19 

m 

■64185 

* Kx) 

2*3 *83291 

1*6 

'47000 

2*0 

'69315 

2*4 *87547 

1-7 

*53063 

2*1 

*74194 


13 

*58779 

2*2 

*78846 
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INTERPOLATION 


2. State which of the two following tables is suitable for linear inter¬ 
polation and estimate the sub tabulation needed in the other to make it 
suitable. 


(i) -176809 

'177840 
■178872 
■179904 
■180937 
■181970 
■183003 
■1S4037 


3, Find /(4*2) from the table. 


X 


X 

2 'G 

•500000 

4-0 

2'5 

■400000 

4*5 

30 

■333333 

5*0 

3'5 

•285714 

5-5 

. The following 

extract is t. 


in Castle *s Five-figure logarithms 


(ii) ■176916 
■177896 
■178891 
479904 
■180937 
481993 
-183075 
484186 


fix) x fix) 

'250000 6-0 466667 

■222222 6-5 453846 

■200000 
481818 

i from the table of natural tangents 
i other tables. 


0 ' 


6 ' 12' IS' 24' 30' 36' 42' 48' 54' 


67 2-35585 36733 37891 39058 40235 41421 42618 43825 45043 46270 

68 247509 48758 50018 51289 52571 53865 56170 56487 57815 59156 


Sf mean differences 

1 2 3 4 3 

67 199 397 596 ^795 994 

68 mean differences cease to be suf¬ 

ficiently accurate 

Find whether linear interpolation does in fact provide 5D accuracy, and 
determine the errors which will be introduced if the mean differences 
are used in calculating the tangents of 67° 3', 67° 27' and 67° 57'. 

5. Evaluate tan \fx when # = *075 radian, using the table for tanV* 
given in problem 4 * 

6 , Subtabulate the following table to produce a table of fix ) for 


*=0(-l)2-5. 

* fix) 

X 

fix) 

X 

fix) 

0 

1-00000 

1*8 

3-10747 

2-3 

503722 

0-9 

1-43309 

1-9 

3-41773 

24 

5-55695 

14 

1-6 

2-15090 

2-57746 

2*1 

2-2 

4*14431 

4-56791 
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EXERCISES 


7. Show that the remainder term for Bessel's interpolation formula 
can be written as 


0(0*_l»)(fla_2 9 ) ... ((fl — n~\ 2 )(8 — n) 
{2 m) I 


(£), 


and find the corresponding remainder terms for the Everett and Stirling 
interpolation formulae. 
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CHAPTER FOUR 


Differentiation and Integration (Quadrature) 


The process of interpolation consists essentially in fitting 
a polynomial to neighbouring tabular values of the func¬ 
tion, and then in substituting the interpolation polynomial 
for the known or unknown function which the table of 
values is supposed to represent. The error of approxima¬ 
tion involved in using P n (x) in place of the actual function 
f{x) can be kept within any prescribed limits over the whole 
of a (finite) arrange* if we use polynomials of sufficiently 
high degree, which in practice implies retaining a sufficient 
number of terms in an interpolation formula. The general 
character of this process is illustrated in figure 2 where the 
continuous curve, which represents an interpolation poly¬ 
nomial Paf#), oscillates about the dashed curve representing 
the true function f(x) and intersects it at tabular values 
(marked by crosses). The lateral separation has been mag¬ 
nified so that the figure may be clearer. If we wish to form 
derivatives or integrals from a table of values we have no 
option but to base our numerical processes on an approxi¬ 
mating polynomial P**(jr) (or some other approximation); 
otherwise we cannot find the derivative, for example, even 
at a tabular point. 

Although the difference between P»(#) and f(x) will be 
small over the working range of the approximation and 
zero at the tabular points, the gradients P^(:c) and f(x) are 
likely to differ more markedly, particularly near tabular 
points. This effect can be illustrated very simply by the 

• Thus the approximation is uniformly valid over the range. 
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non-polynomial (and in fact quite useless) approximation 
S(x) =f(x) e sin(27r xjh) to /(a), where c is small and h is the 
tabular interval. The greatest difference | S(x) —f(x) \ muz 
between f(x) and its approximating function is e at each 
halfway point. However the greatest difference between 
the gradients |3'(*)-/'(*)| max- J(2tt€/A) cos (2frxjh)\ max 
is 2ir«?//t at tabular points, and between Ath derivatives 
/ <fc) (#)l max is {2itjh) k € at tabular points for 
odd-order and halfway points for even-order derivatives. 
Thus for A —-1 the approximate first derivative is in error 
by about 60 c, and the Ath derivative by about (60) I 
Thus numerical differentiation is a process in which 
accuracy is lost > and in particular the determination of 
e 57 

















DIFFERENTIATION AND INTEGRATION 

higher derivatives numerically is usually unsatisfactory 
and is to be avoided if any alternative process is available 
(as by successive differentiation of the function before 
insertion of numerical values, or by successive differentia¬ 
tion of a differential equation, so that higher derivatives 
may be calculated directly by inserting numerical values of 
quite low order derivatives, etc.)* 

The numerical evaluation of an integral J f(x)dx is often 

termed quadrature to distinguish it from the integration of 
differential equations. The integral over any interval can 
be interpreted as the area under the curve in this interval, 
and it follows that integration involves an averaging of the 
errorin each tabular sub-interval and very fre¬ 
quently substantial cancelling of error contributions from 
adjacent sub-intervals (as within the working range of a 
polynomial approximation the error is alternatively +ve 
and ~ve between successive tabular points). When we 
accumulate contributions to an integral over many sub¬ 
intervals the overall accuracy will be greatly enhanced 
because of alternations in sign of local error terms. Thus 
quadrature is a process in which accuracy is gained (in 
powerful contrast with numerical differentiation), and this 
effect is sufficiently marked in practice for us to adopt a 
basic integration process that depends on quadratic (i.e. 
parabola) matching. 

4*1. Numerical differentiation 

Formulae for numerical differentiation can be obtained 
easily by differentiating the interpolation polynomial, and 
this can be carried out conveniently using one of our inter¬ 
polation formulae. Derivatives are needed most commonly 
at a tabular point, and hence from §3.4 we select Stirling’s 
interpolation formula as the form symmetric about tabular 
points. When this is differentiated, 
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4- 1 


303-1 „ . 
+ 3j + 


+ 2M 81/o+ ...}, 

, 1 <PJ\ , 68 ,,, , 2(682-1).,, , 

1 - — 8/ ” + 

+ !2», s%+ )_ 


and so on for higher derivatives, where jc = ,To+0h. These 
forms can be used directly to obtain derivatives at non- 
tabular points by suitable choice of 8 ; and with 0 = 0 we 
obtain finite differmce formulae for derivatives at tabular 
points as, 


hfo'—h^fo- 

- Is4/ 0 +l*y, - ± 8B/ 0 + .... 


It must be emphasized again that these are not infinite 
series, but that the computer must select the truncation 
point appropriate to each application. Thus if a poly¬ 
nomial representation of degree five is appropriate for 
some application we use the form 

which can also be written in the form :* 


/o'^s/o-^sy„+~MS 5 /oj+ 0{h% 

The advantage of the latter form is that it gives an tm- 


• This form implies that in the ease of a function suited to polynomial 
representation the first term neglected (i*e. ^$ 7 / 0 /14OA here) in of order 
h*, and gives a good estimate of the total error for small h. 
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DIFFERENTIATION AND INTEGRATION 

mediate indication of the effect of a change in vvorking 
interval on the truncation error. 

It should be noted that these differentiation formulae 
are of the general form /<*>^ (some combination of dif¬ 
ferences)!}^. The rounding errors in the (combination of 
differences) depend very little on the interval h, so that the 
consequent error in the derivative is proportional to A"*, 
and will be greatly amplified even for low-order derivatives 
when h is smalt. Hence in differentiation we should use as 
large an interval as is consistent with a proper representation 
of the function, even though this may involve the introduc¬ 
tion of a number of higher-order difference terms. 

Example 9. Find d {cos x)jdx at x= 0-5 radians from a 6D 
table of cos x with interval O'OOl. 

From our experience with interpolation we might choose to 
work with as small a tabular interval as possible, -001 in this 
case. Extracting that part of the difference table actually needed, 


X 

/(*) 

¥ 

sy 

*499 

■878062 

-479 

-1 

•600 

*877583 

-480 

-1 

*501 

*877103 


-1. 


The general formula for the first derivative is 

l^s*/o- .... 

and as third differences are negligible we shall truncate it after 
the first term, so that on substitution from the table *001 /'(■ 5) 
— — ’0004795) ^nd f'('5i) — —'480. Note that the maximum 
error in fiSf from rounding of the tabular values f(x) is half a 
unit in the sixth decimal, and that in consequence/'(•&) has a 
possible error of half a unit in the third decimal place 1 

Repeating the calculation with tabular entries at interval -01 
selected from the table, 
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NUMERICAL DIFFERENTIATION 


•49 

•886995 


-88 



-4750 


•50 

■877583 


-88 



-4838 


*51 

•872745 


-88 ; 


we again use the formula A/o— /*8/o and hence -01/'('5) = 
— ‘004794, and /'(■ 5) —— ‘4794 with a possible round-off error 
of half a unit in the fourth decimal digit. Again with A=0*1, 


■4 

■921061 


-9203 


93 



-43478 


434 

-6 

*5 

•877583 


-8769 


88 



-52247 


522 

-4 

•6 

*825336 


-8247 


84 


=—0478626- -0000797- 000000s, and /'(*5)=- *47942 
with a possible error approaching one unit in the fifth decimal 
(arising from the round-off error in the two terms of the formula 
A/o—/di/o-£feS 3 /o). Finally with A=0*3, 


•2 -980067 
*5 -877583 
•8 -696707 



-87547 


7822 

-102484 


9155 



-78392 


7002 

-180876 


16157 



-62235 


5560 


-820 

-1442 


-622 


•3/'(-5)=^-*141680--002109a--0000373+ ... . Still further dif¬ 
ferences must be taken into account (i.e. a higher degree match¬ 


ing is needed), and the difficulty is that at the term ^-^6 5 /o 

uU 


the contribution of the probable rounding error is about 4 x 1Q-? 
and is increasing more rapidly than the coefficients of the series 
decrease; hence we have reached a stage at which the rate of 
decrease of higher-order differences is too slow for numerical 
differentiation (he, the tabular values utilized do not define the 
function sufficiently well for this purpose). 


Note: (i) that the coefficients of the differentiation series 
decrease rather slowly; (ii) that there arc two distinct dis¬ 
advantages of numerical differentiationj the amplification of 
round-off when the tabular interval is small and the succes¬ 
sively poorer matching of the derivatives P^(x) and /**>(#) as 
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the order k increases; (iii) that the calculation of/'(-5) from the 
above data goes much worse, (•001) 2 /''( , 5) ——'000001* 
(■01) 2 /*(-5)^- -000088, etc. 




4.2. The integration of a specified integrand; Simp¬ 
son’s rule 

The process of quadrature, that is the numerical evalua¬ 
tion of the definite integral f{x)dx of some specified 

/unction f(x) t is of considerable importance, as we are often 
unable to find effective analytical processes for evaluating 
integrals even though the functional form of the integrand 
is known. The obvious way out of this difficulty is to ap¬ 
proximate to the known integrand f(x) with some function 
we can integrate, and we are immediately on familiar 


4.2 simpson’s rule 

ground as it was for a very similar kind of approximation 
(to the supposed function underlying a set of tabular values) 
that we have introduced interpolation polynomials, and we 
arc really rather good at integrating polynomials! 

The simplest of interpolation polynomials is the chord 
joining two points of the curve y—f(x) t 

y = P 1 (x)=fo + 6df 0 . 

Hence we have the crude approximation 

f* +a +* ri 

f(x)dx^ j Q Pi(ff)db = A| (f Q +8Afo)d0 t 
J J * 0 Jo 

or 
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which is known as the trapezoidal rule, and consists in re¬ 
placing the integral by the area of the trapezium ABCD, 
as shown in figure 3. A substantial interval a^xKb can 
be subdivided into a number of sub-intervals of length h 
in each of which the simple trapezoidal rale is applied, 
to give 

f ° f(x)dx—h(^fo 4/i 4/2 4 *** -f/tt-i + i/n). 

J a«* o 

The trapezoidal rule will obviously give a poor approxima¬ 
tion unless the sub-interval h is short enough to identify 
each time an almost linear region of the curve y — f(x). 

A much better approximation can be obtained by fitting 
an arc of a parabola to the curve y =f(x) at the three 
tabular points xq , #o4ft, Xq + 2H\ that is by using the inter¬ 
polation polynomial Fa(#) in an interval of length 2/l Then 

j J +2lt f(x)dx 

^jy ah Pi(x)dx = h J *{/ 0 + w/o + d$, 
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Of 

f(x)dx^h(fo + 4/1 + /a). 

This is SijwPjokV and it provides a remarkably simple 
and effective process for numerical integration. We shall 
show that the direct application of Simpson’s rule will in 
many cases provide sufficient accuracy, provided that the 
working interval h is sufficiently small; this will usually 
involve subdivision of the whole interval (a, h) of integra¬ 
tion into an even number of sub-intervals, when 

J * f(x)dx^±h(fo+ 4/1+2/2 + 4/3 + 2/4 + +£.). 

A variety of quadrature formulae can be obtained by 
using interpolation polynomials of successively higher 
degree with matching at a successively larger number of 
tabular points, or by suitable combination of other formulae 
in such a way as partially to cancel errors* However, these 
tend to be very largely of theoretical interest, as they in¬ 
volve subdivision of the range into awkward numbers of 
intervals (which are integer multiples of the degree n of the 
matching polynomial jP n (#)), and also they tend to involve 
a variety of coefficients so that they are less convenient in 
use and more liable to blunders* For example, the three- 
eighths rule is obtained using the matching polynomial 
P%(x) fitted at afo+A, xq + 2 A, xo+$h; and Weddle’s 
rule using /*§(*) fitted at xq, ,,,, *0 + 6 h (see exercise 3 
for these formulae). 

We must have an estimate of the error involved in these 
quadrature processes before we can apply them confidently* 
Suppose that f(x) is continuous in the range (xo — h } Xq + h) r 
and in this range has the Taylor expansion 

where Rz* = x 2k f® k )(£)j(2k)\ for some £=£(#) between 
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4*2 

xo and x. Then on integrating over the symmetrically dis¬ 
posed range (jcq-A, xq + ft), 

fx+h ( 1 L2 M 1 }A 1 fz\ 

■ 2 y° + ®: + 5 + ■ + hJ/jH- 

o 0 

Simpson’s rule can be expressed over the same range in 
terms of the derivatives/^ by expanding/ i=/(xq + ft) and 
f-i = /(xq — h) separately in Taylor series and combining 
the two terminated series, 


f 


*+A 1 

/{*¥* ~ \h(f-y + 4/o +/l) 

x-h -J 


+ 4/0 +/o + hf a + ^ j h% + ... + 1 / (2 **( fz)J 

■=2A{/o+}|yf 0 +j^/o‘ 4 >+ ...j, 

where xq — k < £1 < xo and xq ^ £2 ^ x® + /i* By comparing 
these two forms it may be seen that the correction needed to 
Simpson’s rule to give correct quadrature is 


-s/w-arf**- 


... + ° J 

J m-h 

v 


Rzitdx— 


"3(2 

If f<2X)(x) is bounded in (# 0 —h, xo+h) it is easy to show 
that the remainder term in this correction has the value 

— \{k — l)/i 2 * +1 / (2tJ ( f)/( 2k + 1 )! for some f lying in the range 

(*0 — k, *o + h) but not necessarily equal to any of the 
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previous values. Provided that the interval h is chosen so 
that f(x) is well-represented by a table of values at inter¬ 
val A, the first of these terms will dominate the correction 
and yet will itself be smalh Thus we can integrate a very 
wide range of functions with Simpson's rule provided only 
that vve are prepared to use a suitably chosen working 
interval, and in most cases this will involve subdivision of 
the range of integration; the extent of this subdivision 
will be a matter for judgment based on the character of the 
integrand. We shall also define the error = — {correction). 
The accuracy of formulae based on higher degree inter¬ 
polation polynomials can be estimated in a similar way (see 
exercise 3). Formulae based on polynomials of even degree 
are found in general to be more effective than those of odd 


degree ; thus the corrections expressed as single terms, are 
— fPfZ(£)I 12 for the trapezoidal rule, “A 5 / ( n (£)/90 for 
Simpson's rule, and — 3Ay,^(£)/80 for the three-eighths 
rule (where the ranges of integration have been taken as 
A, 2A, and 3A, respectively, and in each case £ is some point 
within the integration range). Note that the accuracy of 
the three-eighths rule is slightly inferior to that of Simp¬ 
son's rule when the same tabular values are used; it is, 
however, useful to complete an integration over an odd 
number of sub-intervals, when the three-eighths rule can 
be used for either the first three or the last three intervals 
and Simpson's rule for the remainder of the range. 

The above integration formulae use tabular values lying 
within and at the ends of the range of integration, and they 
are sometimes known as closed-type integration formulae. 
We can also develop open-type integration formulae, which 
involve only tabular values within the integration range* 
For example, if we use a parabola fitted to the tabular 
values /i, / 2 , and /a to derive an integration formula over 
the more extensive range (xo t xa+4h) we have* 


• Note that $**0 in the interpolation quadratic corresponds with the 
tabular point x 1 = , e 0 +A. 
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J ^ f(x)dx^hj j/i+uh + $e(e - 

= ^(2/ i -/2 + 2/ 3 ). 

The correction term is + 14/; 5 /< 4 >(£)/45. Closed-type inte¬ 
gration formulae will obviously be more accurate than 
open-type formulae over the same range, as they make use 
of more information; hence closed-type formulae are usually 
used for quadrature, but open-type formulae are useful in 
solving differential equations. 

4.3. Quadrature from a table of values 

Formulae for quadrature from a table of finite differences 
can be obtained by integrating the interpolation formulae. 
Thus from Everett’s formula wc have 

j^f(x)dx -1h jf/o +/i) -1(8% + m + 

++m - ^(8«/ 0 +no )+om 

For integration over an extensive range we can add formu¬ 
lae of this type to span the range (using the results 
S 2 /o + § 2 fi = S/s — §/-;> etc.), 

» +/l +/*+... +f«-l + \fn - 

- Ipy* - % - i*m - 

- - <‘ 8 *> + H 

These formulae represent improved forms of the trapezoidal 
rule with corrections which involve the differences near 
ends of the range of integration; they may also be regarded 
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as giving the error of the trapezoidal rule in terms of finite 
differences of the integrand. 

Similarly, from integration of Stirling’s formula, 

=A{|(/b+4/ 1 +/ 2 )-~syi+^s«/ 1 }+0(n 

giving the corrected form of Simpson’s rule in finite dif¬ 
ferences. 

The choice of truncation point in each of the above 
formulae must be made separately in the context of each 
new problem, but the mode of writing the formulae 
adopted here indicates the dependence on working interval 
of the remainder in case of truncation after sixth or seventh 
differences. In practice it may often be sufficient to work 
with only two or three terms. We often omit remainder 
terms when writing formulae; for in fact we seldom bother 
to evaluate them but rather infer their magnitudes from the 
size and sequential behaviour of the terms which we do 
take into account. Thus we judge directly from the table 
of differences where a formula can be truncated without 
serious truncation error, and in so doing we select the 
degree of the approximating polynomial implicitly. 

The formulae of this section (used in truncated form 
with one or more difference terms) take into account some 
tabular values outside the range of integration as well as 
those within and so correspond with the use of higher 
degree matching polynomials; they are known as partial- 
range integration formulae , and for a specified working 
interval they will be considerably more accurate provided 
that f(x) is well-represented (and this can be judged from 
the size of higher differences). 

Example 10. Evaluate f e*dx from the table 
J o 
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4-3 

x n -A i.n i*r 9 ’ft o- 

e* 1 1-6487 2-7183 4*4817 7*3891 12 1825 

* 30 3*5 4-0 4*5 50 

20*0855 33*1154 54*5982 90-0171 148-4132 

The integrand e % is a smooth function and should be easy to 
approximate except perhaps at the upper end of the range, 
where it is increasing more rapidly. The correct value of the 
integral to 4D is e 5 —1= 147*4132, The trapezoidal rule with 
interval h =0*5 gives 

J e*tf*^-5(|x 1 + 1*6487+2*7183+ 

+ ix 148*4132) 

= 150*4716, with error 3*0584. Simpson’s rule with the same 
working interval gives f e x dx^(-5l'3)[l+ 6*5948+ 5*4366+ 


+ 90*0171+ 


+ 17*9268 + 14*7782 + 48*7300 + 40*1710 + 132*4616 
+ 109*1964 + 360*0684 + 148*4132) = 147*4628, with error 
*0496. The corrected trapezoidal rule needs finite differences 
calculated near the ends of the range, given in the following 


partial table: 






- -5 

■6065 

■1549 


*0392 

*0113 



■3935 


*1003 

-0262 


*0043 

0 

I-0G00 

*2552 


*0654 

■0156 



■64S7 


*1657 

*0413 


■0127 

*5 

L6487 

4209 


■1072 

*0283 


4-5 

900171 

22-9772 


5*8643 

14995 



58*3961 


14-9054 

3*8059 


•9669 

50 

1484132 

37-8826 


9-6702 

2*4664 



96-2787 


24*5756 

6-2723 


1*6022 

5-5 

244-6919 

624582 


15-9425 

4-0686 



Hence 


Pe*<fc — 
J o 


150*4716- (*5/l2){77*3374— *5211) + 


+ (*5 x ll/720)(19*7405 - -1330) - (*5 x 191/60480) x 
x (5*0391-*0340) + (*5 x 2497/3628800X1-2846 - *0085) 

= 147*4132. 

Note: (i) that the simple trapezoidal rule is most seriously 
in error near the upper end of the range where e x varies rapidly 
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with x and there is a contribution 1*9464 from j e x dx to the 

whole error 3*0584; (ii) that even with so large a sub-interval 
as h— 0'5 the parabola fitting of Simpson’s rule gives very rea¬ 
sonable accuracy when the function is smooth in this sub¬ 
interval; (iii) that a standard method of checking an integral 
evaluation is to repeat the calculation using a smaller interval 
(say half of that used initially) and when the two results agree 
to the accuracy required we can be confident both that the 
original interval was suitable and that the calculation is free 
from blunders; (iv) that when greater accuracy is needed the 
finite difference formula should be used; (v) that when the 
behaviour of the integrand varies appreciably from part to part 
of the integration range the most effective procedure may be to 
use working intervals of different length in different parts of the 
range so that the local accuracy of working is roughly uniform; 
(vi) that blunders are specially likely near points at which the 
working interval is changed, because of the interruption of a 
smooth routine, and as a safeguard the two calculations should 
be overlapped to give a check. 


Problem 4, Evaluate 



tan V x dx to four decimals* 


Our first step might be to tabulate the integrand and its 
differences at interval h ={M (see next page). If we evaluate the 
integral by a normal application of Simpson's rule with work¬ 
ing interval /j = 0-1 we have from the table, 
fi .i 

I ^tan Vxdx^-£Q+4x -3272+8 x -4796+ ,..f4x 1*3945+ 

+1-5574)—08536. 

This seems straightforward enough, and as a check on our 
working we repeat the calculation with A=0'05 and we find 
the value 0*8553. This poor agreement is not due to a blunder* 
It can be seen that differences near the top of the table are large 
even at high order, and this is related to the fact that the deriva¬ 
tives of tan V* are all infinite at x^Q. In fact all integrands 
which behave like x l f* or near a terminal point of an 
integration range (or like \x-al±l near x=a) are troublesome 
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X 

V* 

tan y/x 

8 

82 

83 

8 < 

0 

0 

0 

3272 




•1 

•3162 

•3272 

1524 

-1748 

1528 


■2 

•4472 

•4796 

1304 

— 220 

145 

-1383 

•3 

•5477 

•6100 

1229 

-75 

62 

-83 

•4 

•G325 

•7329 

1216 

-13 

38 

-24 

•5 

•7071 

■8545 

1241 

+ 25 

30 

-8 

•6 

•7746 

•9786 

1296 

55 

25 

-5 

•7 

•8367 

1-1082 

1376 

80 

31 

+ 6 

•8 

•8944 

1-2458 

1487 

111 

31 

0 

•9 

•9487 

1-3945 

1629 

142 



1-0 

1-0000 

1-5574 






even though their integrals converge* This kind of difficulty can 
be sidestepped in various ways* 

(a) We can use an open-range formula for one set of intervals 
adjacent to *=0 (allowing us to ‘stand back* a little from the 
trouble level Q) and Simpson’s rule in the remainder of the 
range* Thus with h— 0*1, 

f tan V* d x ~ f tan V x dx+ f tan Vx dx 
Jo Jo J *4 

^=(4 x -1 x JX2 X -3272- -4796+ 2 x -61 G0)+ 

+ (T x i)(*7329+ 4x *8545+ 2 x *9786+ ...+ 

+ 1-5574) 

= 0*8582. 

The correction terms —/ (4) (£) x 10^/90 for each application of 
Simpson’s rule can be seen as negligibly small, but there is 
doubt about the term -f 14/ <4) (£)x 10“ 5 /45 for the initial open- 
type formula. Probably the simplest way of resolving this doubt, 
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and at the same time checking the calculation, is again to repeat 
the process with working interval 0*05, and the "improved 1 result 
is 08569. This decrease in working interval has clearly de¬ 
creased the error from the open-range part of the integration 
but we cannot be sure that it has eliminated it. All that we can 
reasonably assert at this stage is that the integral lies between 
0*8553 and 0*8569, as the correction terms are of opposite sign. 

(b) In the neighbourhood of jc=Q the integrand behaves as 
V*, and hence if we compute the integral of (tan y#—V#) 

and subsequently use the result y#d#=2/3, the numerical 

J o 

work should be a good deal easier. Using values from the table 
above, we have from Simpson's rule with A =0*1, 

I (tan yi- y*)dx=^1895, so that | tan y# dk^'8562. 

Jo Jo 

When the integral is re-evaluated using working interval 
A=G*05 the result is again 0*8562, and this may be taken as 
correct to 4 D. 

(c) An important technique for dealing with awkward 
integrals is to convert them into differential equations in w'hich 
the dependent variable is free from the singularity that has been 

causing the trouble. If we take f tan Vu du=g(x)y(x), where 

J o 


^y(x) is a function to be determined and ^(ar) is selected as having 
a singularity related to that of the integrand, we obtain a dif¬ 
ferential equation on differentiating. With (same 

singularity as tan y# at a;=0) we find 2xy , +y—2V* tan y£; 
with g(x)=^ /2 (when dgjdx has same singularity as integrand) 
we have 2xy'+3y=2tan yjf/yx. For the solution of these 
differential equations see Numerical Solution of Equations . 

(d) A change of variable often serves to reduce a difficult 
integral to more tractable form. In this case, if we take y#= w, 

tan y# dx— 2 u tan u du f 
Jo Jo 

and numerical evaluation is considerably easier. 


Note: (i) that for integration it is normally sufficient to 
work with four decimals to produce a result correct to 4D; 
(ii) that it is worth tabulating the differences of tan y# as a 
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check on the calculated tabular values; (iii) that care is needed 
in working with differences where there is a sign change that 
might be overlooked; (iv) that the variations in fourth differences 
are not so large as to indicate an error in calculation; (v) that 
a difficulty involved in the pursuit of method (a) to smaller 
working interval for the open-range section is that / <4) (£) will 
at the same time increase rapidly; (vi) the methods of (b) and 
(c) are sometimes known as "subtracting out the singularity' and 
"dividing out the singularity’, respectively. 

oe> 

Problem 5. Evaluate ?r 2 to 6D from the series 7r 2 /6="S'(l/ftf 2 ). 

TO-I 

The process of summing terms of a series is related to that of 
integration. In general we need distinct formulae for the two 
types of process, but when the corrected form of the trapezoidal 
rule is expressed in terms of derivatives of the integrand (in 
place of finite differences) we obtain a formula serving both 
processes, and this is the Euler-Maclaurin formula *, 

... +Ul+lfn)-~h>(f:-fo) + 

+ ^:-^)- 3 -okoW 

+ 1209600 A8 ^" > -/<» 7, )+ ••• • 

For the summation of terms of a scries the Euler-Maclaurin 
formula can be rewritten; 

/ 0 +/ 1 + *** +/» 

o* s(c ' /;)+ 
- •)+ 

• For a derivation of the Euler-Maclaurin formula see for example, B* 
T, Whittaker and G. N. Watson, A course of modern analysis (Cambridge 
University Press), chapter 7 P §3; it can also be derived from the trapezoidal 
formula for an extended range by substituting for finite differences m term* 
of derivatives. 
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Hence, when the general term of a series can be written /(m), 
where f(x) is a function possessing the necessary derivatives 
and which can be integrated, we can form sums of terms of the 
series. In the present case f(m)=mr-2; thus /(*)=ar- 2 and 
/(*)(*);=(— !)*(£+1)To sum the series take h= 1 and 
^=1, when 



" +1 ^ + 1 /_J_ + l) + l( _2 \ 

! ar 2\(«+1) 3+ 1 2 J + 12\ (n+1) 3 17 

1 / 24 24 

720 \ (n+l) 5+ l B 


) 


+ 


and if we take the limit of the sum as n-*a>, 


^»i2 - l + 2.1 2+ 6.1 8+ 30.1 s *"* 

These terms decrease very slowly, and we should have to 
accumulate far too many for the required accuracy. This is 
because /(*>(!)=(—I) ! t and it is apparent that a much 
more effective process can be obtained if we sum the first few 
terms directly and apply the Euler-Maclaurin formula to the 
remaining terms. Thus 

yi_i , i . i_i . 1 _ 

xL,nfi 10 2.10 s 6-10 3 30.10 5 42.10 7 
m -10 

= TOOOOOOoH- *0Q500Q0o+ *00016667- *00000033+ *000000Qo 
= *10616634* 

By direct addition 

» 1 

T 6397677a 

i 


and hence ti^Gx+ 64493407=9*869604. The obvious check 
is by repeating the process, taking a different number of initial 
terms for direct addition, and this must be repeated until two 
successive answers agree. 

Note: (i) that this form of the Euler-Maelaurin formula docs 
not imply that it is an infinite series, and for many integrands 
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it will prove to be an asymptotic approximation in which the 
accuracy improves with each successive term only up to a cer¬ 
tain stage and thereafter decreases progressively; (ii) that the 
Euler-Maelaurin formula can be used in its first form for inte¬ 
gration of a specified function provided that derivatives can be 
formed and evaluated with reasonable ease; (in) that we must 
retain eight decimal digits if after adding 9+ 4 terms (each with 
possible round-off of -5 in the last digit) and multiplying by 
6 we are to achieve a result correct to 6D; (iv) that to obtain 
6D accuracy in rfi by direct addition of terms we should have 
to add some twelve million terms of the scries 1 

4.4. Double and repeated integration 

Double (and other multiple) integrals can be evaluated 
by use of the foregoing general techniques applied to a 
repeated integral. Thus 

and if we take F(y) = J f{x, y)dx, we have from Simpson’s 
rule 

4^(-A) + 4F{0) + F(A)j 

= \h[ jj(x, - h)dx +0 )dx+jj(x, h)dx} 

+ 4/o-l +/l-l] + gA[/-10 + 4 /qo +/io] + 

+^(/-u+4/oi+/u]} 

= +/-11 +/-1-1 +/l-l + 4/iO + 4/oi + 4/-10 + 4/0-1 + 
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where In principle this is a perfectly satis¬ 

factory form of Simpson's rule for a double integral, in 
which the value of the integral is given in terms of tabular 
values at the junctions of a square mesh of side-length 2A, 
as indicated in figure 4* However, it will obviously be 
laborious to apply, and we should prefer a formula with 
fewer terms on the right-hand side if this could be obtained 
without undue loss of accuracy. 



One simpler formula, in which the x,y symmetry is 
preserved but from which the relatively distant Corner 
terms 1 /n, /L Uf /_ Wl /i^i are omitted, is of the general 
form 

y) dxd y~ a U +*(/io +/01 +/-10 +/0-1), 

where a and b are constants to be determined. These con- 
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slants can be found easily if we make use of the idea of 
polynomial approximation in a slightly different way, and 
in particular of quadratic approximation which we have 
already found to be effective in integration. If the formula 
that we are constructing is to correspond with quadratic 
matching of the integrand f(x t y) then it will integrate 
correctly (Le. without any error) all quadratic functions 
of x and y f including of course linear functions and con¬ 
stants, Moreover, since integration is a linear process, vve 
can be sure that all quadratic functions will be integrated 
correctly if this is true independently for x 2 , x (and 
hence by symmetry for y 2 t y) t xy and 1, 


When /(#, y) = l we have from the formula 
J a/ ldxdy^xa+4b; 

this double integral can be evaluated directly and has the 
value 4 h 2 t and hence the formula integrates constants 
correctly if a 4 Ah = 4A a . For f(x 9 y) ~ 

0 = J ^xdxdy^axQ + bxQ, 


so that the formula integrates x r and indeed xy and all odd 
powers of x or y (and products odd in x or y) correctly. 
For f(x f y)=x 2 t 

(2A*/3)2A = 

and the formula will integrate *4 correctly (and also y 3 , by 
symmetry) if A = 2A 2 /3, Hence the formula will integrate all 
quadratic functions correctly if a- 4A 2 —4A = 4A 2 /3, and we 
have the double integration formula based on quadratic 
matching, 

a/ ^®( 2 /oo +/i° +/01+/-10 +/0-1) 

**= £(area of element)( 2 /oo +/10 +/oi +/-10 +/0-1), 




ixO+ix 2k 2 , 
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Other formulae of this general type can be obtained in a 
similar way, and these can be tailored to fit any special 
requirements of the problem in hand; indeed, the deriva¬ 
tion above has been given in detail to demonstrate the con¬ 
struction of a formula as well as for its own importance. In 
each case we can calculate a correction term provided that 
a Taylor series (in x and y) is available, and the correction 
terms for the formulae of this section are 

-W+5)tf««(f, V)) 

and 

— 7 })^-fvvvy(^ ^?)}> 

respectively, where (£, is a point of the square domain 

of integration and etc, 

4.5. Gaussian quadrature 

Our previous integration formulae have all been based 
on the use of equally spaced tabular values of the integrand. 
Such formulae are convenient for integration from con¬ 
stant interval tables, and also for integrands of given func¬ 
tional form (as the tabular points can be chosen for easy 
calculation), but they are not the most accurate type of 
formula obtainable. If instead we specify in advance the 
number of tabular values that we wish to incorporate in an 
integration formula but not their position t we can construct 
a formula by finding the tabular positions within the range 
of integration for the most accurate integration formula 
obtainable with the specified number of values of the inte¬ 
grand, This approach is due to Gauss. 

Thus, the most general form of two-point formula of 
this type is 

fj{ X )dx^h{af{cJi) + bm)} t 

where the weighting constants a and b, and the position 
factors a and arc constants to be determined. The approxi¬ 
mation symbol — has been used in anticipation of the fact 
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that we shall now obtain formulae based on polynomial 
approximation (using the device of the previous section). 
The simplest formula of this type* is symmetrical about 
x = 0 and uses values of the integrand calculated at points 
±ah\ for this case /3= We now obtain the values of 
a, b, and a by stipulating that the formula will integrate 
1, x , and x 2 (and hence all quadratic functions) without 
error. 

The constant 1 will be integrated correctly if 2h=h(a + A); 
and x correctly if 0 = h(aoJi — back) ; and x 2 correctly if 
2A 3 /3 =* h(aa 2 h 2 +b<x 2 h 2 ). Hence a = b = 1 , and ofi = 1/3 ; and 
we have the Gauss two-point integration formula 

jj(x)dx^h{f(klV 3) +/(- hj Vi)}; 

this formula (like Simpson’s rule) integrates powers up to 
x 3 (and all odd powers of x) correctly. 

The Gauss two-point formula integrates quadratics 
correctly and corresponds to quadratic representation of 
the integrand /{#), The error may be calculated using 
Taylor expansions, and the leading correction term is found 
to be +(A S / 135)/ (4J (£), which is very similar in magnitude 
to that for Simpson’s rule but is of opposite sign. It should 
be noted that the Gauss two-point formula achieves as 
good accuracy as Simpson’s rule in spite of the fact that it 
employs only two instead of three tabular values; this is 
the advantage of Gauss type formulae, but it must be com¬ 
pared with the disadvantage that these values must be 
determined at awkward points. Thus when we are working 
by hand the inconvenience of evaluating the integrand at 
the Gauss points will more than offset the improved 
accuracy in general and Simpson’s rule (say) will be more 
effective; but when we arc working with a high-speed 
computing machine which can handle any tabular points 


• For special applications we could easily develop non-symmetrical 
formulae of this type. 
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with equal ease the Gauss formulae will have real ad¬ 
vantages. 

The Gauss three-point integration formula 

=l h [ Sf (- h M + 8 ^°> + 5 /( h 7f)) +(/‘715750)/«>(£) 

can be obtained in a similar way, as the formula of type 
h{ a f { *“-f af{pJi) 4- i/(0)} corresponding to cubic approxi¬ 
mation to /(#). 


EXERCISES ON CHAPTER FOUR 

1. Derive the half-way differentiation formula 

s/ . 4 - 2 X + ® 8 \t4 8 \ +0( * , >' 

Find /'(*S5> from the following table for /(x) = <r* and check the result 
by interpolation. 

* -3 -4 -5 -6 -7 8 

/(.v) -74082 -67032 -60653 -54881 -49659 -44933 

2. Obtain the d iflfcrentiation formulae, based on tabular values, 

12A/ 0 '- -/,+8/ 1 -8/_ 1 +/_ 1 +lAyW(f). 

3. (i) Derive the three-eighths rule for integration 

J s h 

Ax)dx^ih(f a + Z/ l + 3/,+f t ) 

using the matching polynomial P,(x), and show that the correction term 
3 

can be written as (£), 

oU 

(ii) Derive Weddle’s rule 

f ^ 3 

J $f(x)dx — Jq^Oo + 5/i H"/i+V 1+/4 + 5j t 
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from the polynomial P/tf). (To obtain this form a small term must be 
neglected-) Find the correction term. 

(iii) Compare the accuracy of Dufton's rule 

5 

I /t^-S:jfc(fi+/*+/i+A) 

frith that of the trapezoidal rule and Simpson’s rule* 


4. Derive the integration formula for use with divided differences 


f f(x)dx^(x t —x 0 )f(x 0 )+fe^o) 

J Xq 

fx i 

show that the correction term is (*-*<,)(*—*i)(*— x^lxx^CiX^dx. 

J 

5. Evaluate f{x)dx from the following table for the error function 


/(*) = V2/7re comparing the results obtained using the trapezoidal 
rule, Simpson’s rule, the three-eighths rule and Weddle’s rule, 


X 

m 

X 

fix) 

X 

/(*> 

►000 

•79788 

*625 

•65632 

1-250 

•36530 

-125 

•79168 

*750 

•60227 

1*375 

•31002 

■250 

•77334 

'875 

■54411 

1*500 

•25903 

■375 

•74371 

loco 

■48394 



■500 

*70413 

1*125 

•42375 




6. Integrate 1/(10-He*+>*) correct to 5D over the square region with 
comers (±1, +!)- [Carry out the integration with the region divided into 
one, four, and then nine elementary squares, and compare the values ob¬ 
tained from these three different working subdivisions of the region.] 


7, Derive the repeated integration formula 
pft 


j f*\ — 2* 1 (Vo+ 3^/0 + j 2 S*/ t jgo/^/o 240^"’" ‘"j’ 

Hence deduce the repeated integration formula 

Jl*/?** - + mro 5 ^ + ■••) * 

8. Construct a double integration formula of type 

^ ^ ^ — a /do +Wu "J"/-n +/- 1 - 1 +/ 1 - 1 )- 

What kind of polynomial matching does it correspond to? Give a Taylor 
series estimate of the correction term. 

81 














CHAPTER FIVE 
Method of Least Squares 


Finite difference processes may be used effectively with 
tables of well-represented functions because the corre¬ 
sponding differences decrease rapidly with increasing 
order. Thus higher differences of mathematical tables 
(with reasonably short interval) tend to be small and regu¬ 
lar in behaviour; and indeed we detect occasional blunders 
by the relatively large disturbances they cause in the other¬ 
wise systematic behaviour of the difference table. The same 
processes will generally be ill-suited for use with tables of 
experimental results, as these will include random fluctua¬ 
tions that often dominate the higher differences. Most 
experimental values will be in error (by an amount which 
exceeds normal round-off), and the effects of these errors 
in the higher differences will be so intertwined that there 
is no hope of detecting (let alone of correcting) them by the 
method of §2.5. We must develop some other way of 
handling experimental results. 

So far we have selected our approximating polynomials 
to agree exactly with the specified tabular values over a 
chosen range of tabular points; thus we tolerate some error 
of approximation between tabular points, though none at 
the points themselves. But when the tabular values include 
random errors there is no advantage to be gained from fit¬ 
ting them precisely, and we would prefer approximating 
polynomials that follow the general trend of a table without 
reproducing all the local fluctuations. This relaxation of the 
accuracy with which we fit individual tabular values in¬ 
volves a reduction in constraint on the approximating poly¬ 
nomial and a consequent increase of range over which a 

82 


FITTING STRAIGHT LINES TO DATA 


S-i 

polynomial of given degree may be used in approximation. 
Matching processes which involve some averaging of small- 
scale local fluctuations are known as curve-fitting, 

5.1. Fitting straight lines to data 

It is often obvious from graphical inspection that the 
measured values, say yt =y(x (), obtained from an experi¬ 
ment lie roughly along a line. These results apparently 
define a linear relationship of the form 

y z±a 0 + aix = Pi(x), 

where the linear representation Pi(x) is not a Newton poly¬ 
nomial (which would have zero error of approximation at 
tabular points) but is some sort of averaged representation; 
we want a procedure for selecting such a ‘line of best fit’ 
or ‘average line’ through the points. It should be noted that 














METHOD OF LEAST SQUARES 

there is no ultimate best line t unless the points are collinear, 
and we shall obtain different ‘best* lines according to the 
criteria of selection we adopt. A good line will pass close to 
most points, and we can select such a line if we adjust 
and ai until the sum of the squares of the residual distances* 
yt-(aq + aixt) is a minimum. The use of (distance) 3 rather 
than simple distance provides a weighting towards distant 
points so that points already near the line exert relatively 
less influence and contributions from points on opposite 
sides of the line reinforce. This criterion of least squares will 
ensure that the distribution of points as a whole lies closely 
about the best line; it is not the only criterion that could 
be used, but it is convenient in many applications. The 
method of least squares can best be illustrated by an example 
of linear approximation. 


Example 11. Use the method of least squares to fit a straight 
line to the following data. 

* 0 1 2 3 4 5 

y 0*82 3*46 6-92 8*57 1107 13-64 

1 he sum of squares S of the residual distances relative to the 
line y=a$+axx is 

5«(0-82—(3*46- oo—ai) 3 +(6-92—*o —2*1)®+ ... 

+ fl3'64—o 1 ®— Siii) 2 * 

In order that 5 may be a minimum for variations of a$ and a x 
(t.e. for variations of the line) it is necessary that 



* — 2{(0'82—flo)+(3-46—<7o—<Ji)+ ... +(13-64-ao-5«i)}, 



= — 2{(3-46-oo-a:)+2(5-92-<io-2ai)+ ... 

^ hence th* +503-64-«,-5« 1 )} : 


• Note that these residua] distances are measured parallel to the y-axii 
or ease, but are proportional to distances from a Line 
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5.2 APPROXIMATION BY LEAST SQUARES 

6i*o4" 43*48 

15a*) 4- 55a* = 153*49* 

Thus the line of best fit in the least squares sense is 
y^G*85+2*56*. 

Note: (i) that we have not proved that this line corresponds 
to a minimum of S t though this is not difficult to do; (ii) that the 
approximate values at tabular points are 0*85, 3*41, 5*97, 8*53, 
11*09, 13*65; (iii) that a measure of the effectiveness of the 
approximation is provided by the root mean square of the 
remaining residuals 

{i[(*03>+(-05)®+(-05)®+ (-04P+ (-02)®+ (-01)®]}*= 0-037; 

(iv) that we cannot strictly determine the accuracy of the 
approximation as we do not know the correct tabular values; 

(v) that we can fit a Newton polynomial P x (x) to just two tabular 
points, but that a least squares polynomial P x (x) can be fitted to 
any number of points. 

5*2. Polynomial approximation by least squares 

In very much the same way we can find polynomial 
approximations /*(*) of any degree for given data; the 
choice of degree is in no way regulated by the number of 
data points, though the accuracy of approximation will 
depend on our choice. The procedure can conveniently be 
described in terms of the algebraically simple case of fitting 
a parabola ^ +atffi to the n +1 tabular values 

f(Xi) (n> 2, and f=0, l t ... ft). 

We must choose the coefficient parameters oq, au and 
az so as to obtain a minimum value for the sum of residualr 

n 

5 = 2 (f(xt)-ao-aiXi - a 2 x$ z . 

t-0 

A set of necessary conditions for 5 to be a minimum is 
that dSjdat- 0 for & = 1, 2; and hence 

n 

-22 m*.) —ao — 1, 2. 














METHOD OF LEAST SQUARES 


n n 

If we take Sk — y and v k = / these equations 
<-0 i-0 

can be written 

So0Q + Si&i + $202 = £*0 
$100 + $201 + $302 = Vi 
S£00 + S^tii + 54^2 = V2j 

and are known as the normal equations of the system. It can 
be shown* that the normal equations always have a unique 
solution for a o> ai f a% and that this solution does correspond 
to a minimum of S . 

The derivation of the normal equations should be laid 
out systematically* as in the following table: 


1 

X 

*2 

3 ? 

X* 

/ 



1 

Xq 

*8 

4 

*0 

fo 

* 0/0 

*o/» 

1 

Xi 


*1 

4 

h 

*l/l 

xffi 

1 

X2 

v2 

*2 

*5 

4 

ft 

* 2/2 

xf % 

1 

X n 



4 

fn 

Xnfn 

X\fn 

*0 

Si 

S2 

*3 

54 

Vo 

Vl 

Vz 


The values xt and /* are given, and may be entered forth¬ 
with; their powers and products are then calculated and 
entered, and the columns are added to obtain s and v sums. 

Example 12. Fit a least squares polynomial to the following 
experimental data (which may be in error to the extent of about 
*02 per value): 

*0 123466 

y — *01 *37 *71 -98 1*28 1*52 1*69 

A substantial saving in labour can be obtained from the 
change of origin £=*—3: in terms of x the tabular points are 

• See Milne, Numerical Calculus (Princeton, 1949), 
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5.2 APPROXIMATION BY LEAST SQUARES 

symmetrically arranged about 5= 0, and if the tabular interval 
is uniform the odd sums are zero, he. fi = 0 f 13 — 0 * •«. * 

We obtain the normal equations for a linear approximation as 
follows: 


Hence 


and 


1 

X 

X* 

y 

*y 

1 

-3 

9 

—01 

•03 

1 

-2 

4 

*37 

--74 

1 

-1 

1 

•71 

-■71 

1 

0 

0 

•98 

0 

1 

1 

1 

1-28 

1-28 

1 

2 

4 

1-52 

8-04 

1 

3 

9 

1-69 

5-07 

7 

0 

28 

6-54 

7-97 


7flo+00i 

= 6-54 



000+2801 

= 7-97, 



>=s= 0934+0*28#*. 



A measure of the effectiveness of this approximation is provided 
by the r.m.s. (root mean square) of the remaining residuals; the 
remaining rms residual is 0*06, and as this is substantially 
larger than the anticipated experimental error it appears that 
this linear approximation is inadequate. 

The normal equations for the quadratic approximation fol¬ 
low as: 


1 

X 

S3 

S3 

P 

y 


x*y 

1 

-3 

9 

-27 

81 

-•01 

■03 

-■09 

1 

-2 

4 

-8 

16 

•37 

-■74 

1-48 

1 

-1 

1 

-1 

1 

■71 

-•71 

•71 

1 

0 

0 

0 

0 

■98 

0 

0 

1 

1 

1 

1 

1 

1-28 

1-28 

1-28 

1 

2 

4 

8 

16 

1*52 

3-04 

608 

1 

3 

9 

27 

81 

1-69 

5*07 

15*21 

7 

0 

28 

0 

196 

6-54 

7*97 

24-67 


7tf0+ 00i+ 2802= 6*64 
000+2801+002=7*97 
2800+001+19602—24*67- 
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METHOD OF LEAST SQUARES 

Hence y 1-005+0-28s*— 0-0 la* 8 

= -0-01z+0-39s*—O-Ols* 8 . 

The rms residual error is now 0-01, which is well within the 
range of possible experimental error and must be regarded as 
satisfactory. 

Note; (i) that we now rely on our judgment in choosing the 
degree of approximating polynomial that is appropriate, and 
there is nothing to prevent us from finding a ‘best straight line’ 
through points obviously ill-suited for linear approximation; 
(ii) that the root mean square residual is a valuable way of 
assessing an approximation by least squares; {iii) that there is 
an effective reduction by one in the number of normal equations 
to be solved when the tabular points are arranged symmetric¬ 
ally; (iv) that the third decimal digits are needed as subscripts 
(but cannot genuinely be written as ‘correct digits') in order to 
preserve accuracy. 

Experimental results are often plotted on log-log or 
log-plain graph paper to bring out relationships of power 
law and exponential type. The same procedure can be 
carried out much more accurately using least-square linear 
approximation. In the case y — ax n , we have log y — log a 
+ « log x, and if we put X= log x, 7=log y, and A = log a, 
we have, 

Y=A + nX; 

A and « can be determined now by a least-squares process 
applied to the values logjy* at points log x\. Similarly, if 
y = ae bx , we take Y= log y and A = log a to give Y—A + bx, 
and a least-squares process can be applied to the values 
log yt at points xj. 

It should be emphasized that the least-squares approach 
is bound to yield an approximation,* though we shall get 
a useless approximation if we select the degree of the ap¬ 
proximating polynomial foolishly. It follows that we must 
always test the significance of an approximation by a final 
calculation of the rms residual, or in some other way. 

• Compare the comments on 6emi-empiricai theories of §3,7 (iii). 
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5-2 APPROXIMATION BY LEAST SQUARES 

In some applications the tabular values may be regarded 
as of unequal importance or weight-, in such cases a weight 
function w(* ( ) can be associated with each tabular point, and 
the least-squares process will then be carried out on the 
weighted sum 

^w{x t ){ffx t )-P n {xi)}\ 

Problem 6* Find a least squares quadratic approximation for 
the continuous function e x over the range (— 1, 1), 

We might select an arbitrary set of tabular points spanning 
the range (-1, I) and apply the methods of this section to the 
corresponding set of values for But it will be obviously 
superior to use the known values of e x right across the range by 
using integration in place of the summation process necessary 
for sets of discrete values. Thus, in place of the sum of squares 
of the residuals we take the integral 

S= J {e*— (oq + a&+ a 2 x 2 )} 2 dx, 

which is a function of oq, and a 2 - A necessary set of conditions 
for a minimum of the function S(ao t a x , a%) is that dS/da^O 
for k= 0,1, 2: hence 

2 O i*l 

0= — 2J {e x — (flo+ «i*+ 02 -x z )}dx, 

bs n 

0 = = — 2 J {e z — (oq+ a x x + tfs# 2 )}# dx t 

bs r l 

0 = — 2 j (# 0 + ffi *+ Q 2 jc 2 )}* 5 dx. 

After integration, we obtain the set of normal equations 
Stfo-Ma —3*53 
01=1*10 
5<to+ 3^2— 6*59* 

and the quadratic approximation 

y^ 1*00+ HC*+ 0-53*®. 

Note: (i) that the numerical value for e has been substituted 
immediately as we are concerned with numerical solutions; 
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METHOD OF LEAST SQUARES 

(ii) that the error does not exceed about '03 units except at the 
very ends of the range, and that the rms residual error is *03; 

(iii) that the approximation differs from the first three terms of 

the Taylor expansion near x= 0, but is more accurate over the 
specified range £5 (iv) that this process of matching a 

continuously defined function is known as continuous approxi¬ 
mation by least squares. 

5.3. Data smoothing 

Raw experimental results provide our only source of 
contact with the physical world; we should on.no account 
‘adjust’ these results capriciously, either to improve an 
‘apparent pattern of dependence’ (which may not be a 
genuine one, in any case) or to provide a better level of 
agreement with a favoured theory. Experiments which 
yield consistently poor results may need refining, or they 
may be based on an inadequate conception of the pheno¬ 
mena involved; in either case we can improve the apparent 
quality of the results by smoothing processes , but we will not 
improve their significance (though we may appear to do sol). 
Thus data smoothing is not an alternative to good experi¬ 
mental technique. However, it must be appreciated that 
even good experiments yield tables of results with enough 
experimental scatter to produce sizable spurious fluctu¬ 
ations in the higher differences. Some initial smoothing 
may be necessary before such tables are used in numerical 
calculations, especially for processes like differentiation 
which are very sensitive to local error. 

In smoothing a table of values we do not attempt to find 
a single approximating polynomial for the whole table; 
instead we use repeated local matching with low-degree 
polynomials, and each polynomial is used only to adjust a 
single central value (except at the ends of the range). Thus 
we assume that the table is locally suitable for approxima¬ 
tion by polynomials of low degree, apart from the experi¬ 
mental scatter, but we make no assumption about large- 
scale approximation of the table. 


5.3 DATA SMOOTHING 

A wide range of formulae is available for smoothing data 
at uniform interval; we shall derive a quadratic five-point 
formula. We shall suppose that a run of five values is suf¬ 
ficient to correct for random errors in the measured value 
of the central member, and we shall take the tabular 
interval as the unit of length and the central point as 
origin. The normal equations for fitting a parabola 
ao+aix+axt? to the values/-s,/-i, /o, ft, /s by least-squares 
are then 

a 

5ao+10«a = £/i 

—2 

2 

-a 

a 

Hko + 34fl2 = 5>|f*i. 

—2 

For use in the interior of the table we want only a corrected 
value for the central point, that is the value a 0 ; and we 
have on solution 

35do = 35/o - 3(/-2 - 4/-i + 6/0 - 4/i +/a) 

= 35/o —3S 4 /o. 

Hence the corrected value* at a: = 0 is yo — /o — 38yo/35. 
For use near the head or foot of a table 

y-i —f-i+ ^S 4 /o, yi —ft+ 35^/0. 

In selecting the degree n for the smoothing process we 
should take n small so that individual calculations are easy, 
and yet we must so judge our choice that differences of the 

* For many purposes it will be sufficient to use the approximate form 

y # Sy„/12. 
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true function of order greater than n are small. In fact the 
choice may often be made from the fact that the first, n 
differences are reasonably regular, but the (n+l)th dif¬ 
ferences are irregular and have mean value near zero. The 
amount of smoothing increases with the number of values 
fitted, and decreases with the degree of the approximating 
polynomial. Obviously there is great scope lot judgment on 
the part of the computer as to what is random and what 
genuine; sometimes the application of successive smooth¬ 
ing processes may be justified, but repeated smoothing will 
ultimately reduce any set of data to a straight line. 


EXERCISES ON CHAPTER FIVE 


1. Use the method of least squares to find the best straight line fitting 
the data 

* 0 1 2 3 4 56789 10 

y 15 14 13 13 12 10 6 8 5 4 2; 

2. Find the best curve of type y*=ax n fitting the data 

tt 0-13 0-33 1-08 3-65 4-98 
y 3 4 6 9 10. 
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Note on Operational Methods 


Most readers will have some familiarity with operational 
methods through their use in solving linear differential equa¬ 
tions with constant coefficients. It is possible in a purely formal 
way to dissociate the operation of differentiation (or integration) 
from any particular function f(x) and to establish a formal 
algebra for the operator D=dfdx t which is similar to (though not 
identical with) the algebra of real numbers. 

In precisely the same way we can regard the operations of 
differencing as distinct from their application to a particular 
function, in which case A f p and 8 may be regarded as differenc¬ 
ing operators and p as an averaging operator . A formal algebra 
can be established in terms of these and other operators, and 
this operational algebra can be used in deriving a very wide 
range of formulae for numerical processes. Operational methods 
are used in many books because they allow of such easy develop¬ 
ment of formulae, but they have not been used here as they 
tend to leave the significance of the various steps obscure to the 
inexperienced reader. We have advanced more slowly, in the 
hope that the reader will carry a clear understanding of the 
reasons behind our development of processes and of our atti¬ 
tude towards numerical work; he can then take up the opera¬ 
tional approach with help from a more advanced text. 
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Answers and Notes on Exercises 


Chapter One 

1. (i) 3, 6; about 1 b 2x10®, (ii) 30, 3; about Q OS per cent, (iii) 0, 5; 
about 1 in 60,000, (iv) 4, 3; about 0 4 per cent, and note that the final 
zero is needed to define this accuracy. 

2. 1*38 and 1-37; note that the rounding error can pV st ) ex ? c ed half a 
unit in the last digit when a number is rounded digit by digit. 

3. (i) 0-5; actually 0-47+7 X-DOS, and hence the maximum error in the 
second decimal from round-off is +3^, so that this digit is uncertain 
(the probable error is considerably smaller, but exceeds *005). 
(ii) 83-5; maximum error is about (9*374-8-91)X*005—0-09, and the 
second decimal in 83*49 is unreliable and the error could affect the 
first decimal, (iii) 25-6; maximum rounding error is 0 048 after 
division and the second and third decimals are unreliable (clearly^the 
answer cannot be as accurate as the denominator 0*279!). (iv) 17-760; 
the decimal accuracy of these numbers varies and as the middle num¬ 
ber is known only w ithin *0005 there is no justification for retaining 
more than three decimals, maximum error is about *0006. Note that 
maximum errors can be found using the binomial expansion, as 

=|(l±|) (l±|) 1 =|( 1 ±( fl +|)+ -)• 

4. (i) 0*693144 with maximum rounding error of 1 in last digit, (ii) 
0 1 693143 with maximum rounding error of 4 in last digit; this illus¬ 
trates the value of carrying an extra working digit for some calculations* 

5. (i) 14-tf, for the range -*!<*<*!; (ii) + for the range 
—3<ac<-3; (iii) 14 A-4-i* a 4-£a\ for the range — 5<a<*5, 


Chapter Two 

2. Compare w ith formula for differentiation of products, etc. 

3. Use induction. 

4. Third differences are constant: hence data fitted by a cubic which is 
obtained from an interpolation formula as /5 -r521.V; 6-^25x a i2 4* 2.V 3 3* 

5. /(1*S) = 1 142892. Transposition is a common blunder; the adjusted 
fourth differences am about as regular as can be expected, allowing 
for round-off. 
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ANSWERS AND NOTES ON EXERCISES 


Chapter Three 

1. Use Ncw+on-Gregory interpolation with forward differcncea near 
head of table and backward differences near foot. Differencing of final 
table will uncover any errors, so that there is no need to give solution 
here. 

2. (i) Suitable for linear interpolation as S 2 / < 1 over whole table, (ii) Not 
suitable* as 8 s / rises from 15 to 29. Subtabulation to half-interval will 
reduce second differences by factor of about 1/2®; they are still too 
large to be neglected. Subtabulation to one-third interval will reduce 
second differences by factor of about 1/3®, when they may be neglected 
over whole range (as being ^ 4 units in last decimal), 

3. /(4*2) =0238095. Throwback of fourth on to second differences is 
not good enough because of the large interval and number of decimal 
digits; use Everett to sixth differences. 

4. §*/ is about 10, and makes a contribution to Bessel's formula of about 
one unit in the last decimal over much of the 0-range 1 Values for the 
tangents using Besselformula are 2*36158, 2-40827, 2*46888, and 
using mean differences are 2-36181, 2-40831, 2*46866, respectively, 
and thus errors of more than twenty units in the last decimal are possible 
near 0' and 60' 1 In fact the values using mean differences for 67° 3' 
and 67° 57' are correct to 3D only HI and that for 67 5 27' to 4D! It 
is hard to see why the author retains mean differences for 67° but 
rejects them for 68 3 ; and the reader may find it instructive to investi¬ 
gate other short sets of tables and to compare them with Chambers's 
Shorter Six-figure Mathematical Tables, 

5. First reaction is to use the Newton-Grcgory forward difference inter¬ 
polation formula as w r e are working near the head of a table, but we 
find that the terms decrease too slowly for effective calculation. 
{N<B. if we take many terms this will involve many tabular values, 
only a few* of wdnch will be from the neighbourhood of x = 0, where 
tan V*** V*.) We can avoJd ^ is difficulty by squaring each term in the 
table and working with tan 2 \fx (which behaves like x near A?—0>; 
after the interpolation has been carried out successfully for tan a (*075) 
take the square root. 

6. Use divided differences for whole subtabulation, or alternatively to 
tabulate at #=0(-3)2-4 follow-ed by Newton-Grcgory at constant 
interval. Former method probably simpler once divided differences 
have been calculated. Check final table by differencing. Actually 

and solution may be checked from tables, though it is 
essential that the computer should gain confidence in bis own checks 
for ensuring accuracy- 

7. Everett remainder same as Bessel; Stirling remainder term is 

ftp-i*)# 1 -?)... (gw? Srt+1 , (i „ (A 
(2»+l)l 
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ANSWERS AND NOTES ON EXERCISES 


Chapter Four 

1. Differentiate Bessel's interpolation formula; use result /'=* —/ IQ 
check result by interpolation on table given. 


2. Obtained from finite difference formula by breaking down the first, 
and the first two terms, respectively, into tabular values. 

3. (ii) —V/W(f)/140; (iii) — 5/i J /"(£)/6 compared with — Nf'X$)ll2 
for each application of the trapezoidal rule and —A^ 4 l(f)/90 for each 
application of Simpson’s rule. Dufton’s rule is comparable with the 
trapezoidal rule in accuracy and is remarkably easy to apply: it is use 
ful for a quick estimation of an integral from a table. 


5* Correct value is 0*86638, 


6. Correct value is 0*37558, 

7. Obtained by integration of Stirling’s interpolation formula. Second 
form for integration over a double interval obtained by adding 
formulae of first type for 

and last term obtained by suitable modification of original formula* 
M Ah 

8* I dx /(x, y)dy ^ {8/^ -f/ u -f/-n +/- 1-1 +/ 1 - 1 } - 

J -h J -h 

Matching is quadratic, and correction term is 
—45^* {fxxxx (fc*)+l Q fxxyyi f 


Chapter Five 

1. 16— \ '3*. Note that the rms residual error is about 1, and there ia 
nothing to be gained from specifying the line more accurately» 

2. S-SSa 3 * 3 , with rms residual error of about 0 03, 
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Directory of Numerical Processes 


Polynomial approximation 

Exact matching of equidistant tabular values (see Interpolation) 
Least squares matching of discrete tabular values: §§5.1,5,2, 
Continuous approximation by least squares: problem 6 (§ 5.2), 


Other approximations 

Exponential and power law approximation by least squares: 
§5,2, 


Interpolation 

Using divided differences based on table with variable interval: 

§3-2. 

Linear interpolation: § 3,1, 

Using differences of a constant interval table: §§3i4> 3,5 j 
including throwback: § 3.5. 


Differentiation 

Using differences of a constant interval table : § 4.1. 

Using tabular values only: exercise 2 (Ch. 4). 

Half-way formula using finite differences: exercise 1 (Ch. 4). 

Integra tion (Quadra ture ) 

Using ordinates at constant interval without differences: § 4.2. 
Using constant interval differences: §4.3. 

When the integrand has a singularity: problem 4 (§4.3). 
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DIRECTORY OF NUMERICAL PROCESSES 

Double and repeated integration : § 4.4, exercise 7 (Ch* 4), 
Using ordinates at ‘optimum* points: §4.5, 

Using divided differences : exercise 4 (Ch, 4), 

Other processes 

Detection of errors in a table: example 2 (§ 2,5), 

Extraction of roots of an algebraic equation: example 8 (§ 3.6)* 
Interpretation of a number: § 1.3. 

Line of‘best fit* for (experimental) data: § 5*1* 

Smoothing of irregular data: § 5.3* 

Summation of series: problem 5 (§ 4.3)* 
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Index 


Accuracy of a number, 7 

Backward differences, 21 
Bessel coefficients, 42, 52 
Bessel interpolation formula, 
42 , 44 

remainder, 55 
Blunders, 10, 11 

Central differences, 21 
Continuous approximation, 

90 

Continuous functions, 8 
Critical tables of coefficients, 
52 

Curve fitting, 9, 83 
D, 15 

Data smoothing, 90 
Decimal digits, 7 
8,22 
d, 20 

Detection of errors in tables, 
25,26 

Differences, 17, 18 
backward, 21 
central, 21 
divided, 28, 29, 30 

expected fluctuations, 20 
forward, 20 


modified for throwback, 46, 
51 

of a polynomial, 22, 23, 29 
Differentiation, 56, 57 
numerical processes, 59, 80 
errors, 60 

Divided differences, 28 
Double integration, 75, 77 
Simpson’s rule, 75 
Dufton’s rule, 81 

6, 37 

Error fan, 25 
Errors, 10 
detection, 25 

in numerical differentiation, 
57, 60 

in numerical integration, 64 
probable, 10 
random, 9 
round-off, 10, 19 
systematic, 9 
truncation, 11 

Euler-Maclaurin summation 
formula, 73 

Everett interpolation formula, 
42,44 

Exponential approximation, 88 
Extrapolation (Newton-Greg¬ 
ory formula), 40 
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INDEX 
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Gaussian quadrature, 78, 79 

Integration, 58, 62 
divided difference formula, 
81 

Dufton’s rule, 81 
errors, 64, 66 
Gaussian, 78, 79 
Simpson's rule, 64, 68 
three-eighths rule, 64, 80 
trapezoidal rule, 63, 67 
Weddle’s rule, 80 
Interpolation, 8, 32, 33, 56 
Bessel formula, 42, 44 
choice of process, 44,47 
divided differences, 34, 36 
Everett for mu la, 42, 44 
half-way formula, 43 
inverse, 49 
linear, 32 

Lagrange formula, 38 
N e wton - Gr e go ry formula, 
39 

Stirling formula, 43 
Inverse interpolation, 49 
Iterath r e process, 51 

Lagrange interpolation formu¬ 
la, 38 

Least squares, 84 


continuous approximation, 
90 

exponential approximation, 

88 

linear approximation, 84, 

87 

polynomial approximation, 
85, 87 

power law approximation, 

88 

rjn.s* residual, 85 
Linear interpolation, 32 

*42 

Ft 21 

Newton-Gregory interpolation 
formula, 39 

Newton interpolation formu¬ 
la, 36 

Normal equations, 86 
Notation for finite differences, 
20 

for tables, 15 
Numbers, 7 

Operational methods, 93 

Periodic functions, 16 
Pivotal values, 48 
Polynomial representation, 
3-6, 24 

Power law approximation, 88 
Proportional parts, 34 

Quadrature (see integration) 
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Smoothing, 9, 90 
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Summation of series, 73 
Systematic error, 9 

Tables of values, 15 
Tabular interval, 15, 20 
Taylor series, 2, 53 
Throwback, 45 
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Truncation error, 3 

Weierstrass’s theorem, 4 
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